From 4a3d668c59eefa871c668e58afd98a099affdd43 Mon Sep 17 00:00:00 2001
From: Tobias Leibner <tobias.leibner@googlemail.com>
Date: Wed, 26 Jun 2019 09:09:54 +0200
Subject: [PATCH] [test] Fix tests for kinetic scheme
---
dune/gdt/momentmodels/density_evaluations.hh | 172 +++++
.../entropyflux_implementations.hh | 627 +++++++++++++-----
.../momentmodels/entropyflux_kineticcoords.hh | 84 ++-
dune/gdt/momentmodels/hessianinverter.hh | 49 +-
.../operators/advection-fv-entropybased.hh | 11 +-
dune/gdt/operators/reconstruction/internal.hh | 2 -
.../reconstruction/linear_kinetic.hh | 31 +-
dune/gdt/operators/reconstruction/slopes.hh | 26 +-
.../test/entropic-coords-mn-discretization.hh | 78 +--
...bolic__momentmodels__entropic_coords_mn.cc | 58 +-
.../momentmodels/kinetictransport/base.hh | 30 +
.../kinetictransport/sourcebeam.hh | 128 +++-
.../kinetictransport/testcases.hh | 159 ++++-
13 files changed, 1081 insertions(+), 374 deletions(-)
create mode 100644 dune/gdt/momentmodels/density_evaluations.hh
diff --git a/dune/gdt/momentmodels/density_evaluations.hh b/dune/gdt/momentmodels/density_evaluations.hh
new file mode 100644
index 000000000..2e6a212d6
--- /dev/null
+++ b/dune/gdt/momentmodels/density_evaluations.hh
@@ -0,0 +1,172 @@
+// 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_DENSITYEVALUATOR_HH
+#define DUNE_GDT_MOMENTMODELS_DENSITYEVALUATOR_HH
+
+#include <string>
+#include <functional>
+
+#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, SlopeType slope>
+class LocalDensityEvaluator : public XT::Grid::ElementFunctor<typename SpaceType::GridViewType>
+{
+ using BaseType = XT::Grid::ElementFunctor<typename SpaceType::GridViewType>;
+
+public:
+ using GridViewType = typename SpaceType::GridViewType;
+ using EntityType = typename GridViewType::template Codim<0>::Entity;
+ using IndexSetType = typename GridViewType::IndexSet;
+ using EntropyFluxType = EntropyBasedFluxEntropyCoordsFunction<GridViewType, MomentBasis, slope>;
+ using RangeFieldType = typename EntropyFluxType::RangeFieldType;
+ 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>;
+ using DomainType = FieldVector<RangeFieldType, dimFlux>;
+ using BoundaryDistributionType = std::function<std::function<RangeFieldType(const DomainType&)>(const DomainType&)>;
+
+ explicit LocalDensityEvaluator(const SpaceType& space,
+ const VectorType& alpha_dofs,
+ EntropyFluxType& analytical_flux,
+ const BoundaryDistributionType& boundary_distribution,
+ const RangeFieldType min_acceptable_density,
+ const XT::Common::Parameter& param)
+ : space_(space)
+ , alpha_(space_, alpha_dofs, "alpha")
+ , local_alpha_(alpha_.local_discrete_function())
+ , analytical_flux_(analytical_flux)
+ , boundary_distribution_(boundary_distribution)
+ , min_acceptable_density_(min_acceptable_density)
+ , param_(param)
+ , index_set_(space_.grid_view().indexSet())
+ {}
+
+ explicit LocalDensityEvaluator(LocalDensityEvaluator& other)
+ : BaseType(other)
+ , space_(other.space_)
+ , alpha_(space_, other.alpha_.dofs().vector(), "source")
+ , local_alpha_(alpha_.local_discrete_function())
+ , analytical_flux_(other.analytical_flux_)
+ , boundary_distribution_(other.boundary_distribution_)
+ , min_acceptable_density_(other.min_acceptable_density_)
+ , param_(other.param_)
+ , index_set_(space_.grid_view().indexSet())
+ {}
+
+ virtual XT::Grid::ElementFunctor<GridViewType>* copy() override final
+ {
+ return new LocalDensityEvaluator(*this);
+ }
+
+ void apply_local(const EntityType& entity) override final
+ {
+ local_alpha_->bind(entity);
+ const auto entity_index = index_set_.index(entity);
+ XT::Common::FieldVector<RangeFieldType, dimRange> alpha;
+ for (size_t ii = 0; ii < dimRange; ++ii)
+ alpha[ii] = local_alpha_->dofs().get_entry(ii);
+ analytical_flux_.store_density_evaluations(entity_index, alpha);
+ for (auto&& intersection : Dune::intersections(space_.grid_view(), entity))
+ if (intersection.boundary())
+ analytical_flux_.store_boundary_evaluations(
+ boundary_distribution_(intersection.geometry().center()), entity_index, intersection.indexInInside());
+ } // void apply_local(...)
+
+private:
+ const SpaceType& space_;
+ const ConstDiscreteFunctionType alpha_;
+ std::unique_ptr<typename ConstDiscreteFunctionType::ConstLocalDiscreteFunctionType> local_alpha_;
+ EntropyFluxType& analytical_flux_;
+ const BoundaryDistributionType& boundary_distribution_;
+ const RangeFieldType min_acceptable_density_;
+ const XT::Common::Parameter& param_;
+ const typename SpaceType::GridViewType::IndexSet& index_set_;
+}; // class LocalDensityEvaluator<...>
+
+template <class MomentBasisImp,
+ class SpaceImp,
+ SlopeType slope,
+ class MatrixType = typename XT::LA::Container<typename MomentBasisImp::RangeFieldType>::MatrixType>
+class DensityEvaluator
+ : 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, slope>;
+ using RangeFieldType = typename MomentBasis::RangeFieldType;
+ using LocalDensityEvaluatorType = LocalDensityEvaluator<SpaceType, VectorType, MomentBasis, slope>;
+ using BoundaryDistributionType = typename LocalDensityEvaluatorType::BoundaryDistributionType;
+
+ DensityEvaluator(EntropyFluxType& analytical_flux,
+ const SpaceType& space,
+ const BoundaryDistributionType& boundary_distribution,
+ const RangeFieldType min_acceptable_density)
+ : analytical_flux_(analytical_flux)
+ , space_(space)
+ , boundary_distribution_(boundary_distribution)
+ , min_acceptable_density_(min_acceptable_density)
+ {}
+
+ 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& alpha, VectorType& /*range*/, const XT::Common::Parameter& param = {}) const override final
+ {
+ analytical_flux_.density_evaluations().resize(space_.grid_view().size(0));
+ analytical_flux_.boundary_density_evaluations().resize(space_.grid_view().size(0));
+ LocalDensityEvaluatorType local_density_evaluator(
+ space_, alpha, analytical_flux_, boundary_distribution_, min_acceptable_density_, param);
+ auto walker = XT::Grid::Walker<typename SpaceType::GridViewType>(space_.grid_view());
+ walker.append(local_density_evaluator);
+ walker.walk(true);
+ } // void apply(...)
+
+private:
+ EntropyFluxType& analytical_flux_;
+ const SpaceType& space_;
+ const BoundaryDistributionType& boundary_distribution_;
+ const RangeFieldType min_acceptable_density_;
+}; // class DensityEvaluator<...>
+
+
+} // namespace GDT
+} // namespace Dune
+
+#endif // DUNE_GDT_MOMENTMODELS_DENSITYEVALUATOR_HH
diff --git a/dune/gdt/momentmodels/entropyflux_implementations.hh b/dune/gdt/momentmodels/entropyflux_implementations.hh
index 5d3fdbc30..b3f3e78eb 100644
--- a/dune/gdt/momentmodels/entropyflux_implementations.hh
+++ b/dune/gdt/momentmodels/entropyflux_implementations.hh
@@ -50,6 +50,14 @@ namespace Dune {
namespace GDT {
+template <class T>
+std::enable_if_t<std::is_arithmetic<T>::value, T> superbee(const T first_slope, const T second_slope)
+{
+ return XT::Common::maxmod(XT::Common::minmod(first_slope, 2. * second_slope),
+ XT::Common::minmod(2. * first_slope, second_slope));
+}
+
+
// choose specializations
#ifndef ENTROPY_FLUX_UNSPECIALIZED_USE_ADAPTIVE_CHANGE_OF_BASIS
# define ENTROPY_FLUX_UNSPECIALIZED_USE_ADAPTIVE_CHANGE_OF_BASIS 1
@@ -72,6 +80,15 @@ namespace GDT {
#endif
+enum class SlopeType
+{
+ minmod,
+ superbee,
+ mc,
+ no_slope
+};
+
+
template <class MomentBasisImp>
class EntropyBasedFluxImplementationUnspecializedBase
: public XT::Functions::FunctionInterface<MomentBasisImp::dimRange,
@@ -105,7 +122,8 @@ 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>;
+ using QuadraturePointsType = std::vector<BasisDomainType, boost::alignment::aligned_allocator<BasisDomainType, 64>>;
+ using QuadratureWeightsType = std::vector<RangeFieldType, boost::alignment::aligned_allocator<BasisDomainType, 64>>;
explicit EntropyBasedFluxImplementationUnspecializedBase(const MomentBasis& basis_functions,
const RangeFieldType tau,
@@ -230,7 +248,7 @@ public:
const size_t dd) const
{
- thread_local FieldVector<std::vector<RangeFieldType>, 2> work_vecs;
+ thread_local FieldVector<QuadratureWeightsType, 2> work_vecs;
work_vecs[0].resize(quad_points_.size());
work_vecs[1].resize(quad_points_.size());
calculate_scalar_products(alpha_i, M_, work_vecs[0]);
@@ -260,8 +278,8 @@ public:
// 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
- template <typename T, typename Compare>
- static std::vector<std::size_t> get_sort_permutation(const std::vector<T>& vec, const Compare& compare)
+ template <class T, class Alloc, class Compare>
+ static std::vector<std::size_t> get_sort_permutation(const std::vector<T, Alloc>& vec, const Compare& compare)
{
std::vector<std::size_t> p(vec.size());
std::iota(p.begin(), p.end(), 0);
@@ -269,8 +287,8 @@ public:
return p;
}
- template <typename T>
- static void apply_permutation_in_place(std::vector<T>& vec, const std::vector<std::size_t>& p)
+ template <class T, class Alloc>
+ static void apply_permutation_in_place(std::vector<T, Alloc>& vec, const std::vector<std::size_t>& p)
{
std::vector<bool> done(vec.size());
for (std::size_t i = 0; i < vec.size(); ++i) {
@@ -290,8 +308,7 @@ public:
}
// Joins duplicate quadpoints, vectors have to be sorted!
- static void join_duplicate_quadpoints(std::vector<BasisDomainType>& quad_points,
- std::vector<RangeFieldType>& quad_weights)
+ static void join_duplicate_quadpoints(QuadraturePointsType& quad_points, QuadratureWeightsType& quad_weights)
{
// Index of first quad_point of several quad_points with the same position
size_t curr_index = 0;
@@ -319,7 +336,7 @@ public:
static_assert(std::is_same<BasisFuncImp, MomentBasis>::value, "BasisFuncImp has to be MomentBasis!");
RealizabilityHelper(const MomentBasis& basis_functions,
- const std::vector<BasisDomainType>& quad_points,
+ const QuadraturePointsType& quad_points,
XT::Common::PerThreadValue<std::unique_ptr<ClpSimplex>>& lp)
: basis_functions_(basis_functions)
, quad_points_(quad_points)
@@ -386,7 +403,7 @@ public:
private:
const MomentBasis& basis_functions_;
- const std::vector<BasisDomainType>& quad_points_;
+ const QuadraturePointsType& quad_points_;
XT::Common::PerThreadValue<std::unique_ptr<ClpSimplex>>& lp_;
}; // struct RealizabilityHelper<...>
#else // HAVE_CLP
@@ -432,16 +449,16 @@ public:
}; // struct RealizabilityHelper<Hatfunctions, ...>
// temporary vectors to store inner products and exponentials
- std::vector<RangeFieldType>& working_storage() const
+ QuadratureWeightsType& working_storage() const
{
- thread_local std::vector<RangeFieldType> work_vec;
+ thread_local QuadratureWeightsType work_vec;
work_vec.resize(quad_points_.size());
return work_vec;
}
void calculate_scalar_products(const DomainType& beta_in,
const BasisValuesMatrixType& M,
- std::vector<RangeFieldType>& scalar_products) const
+ QuadratureWeightsType& scalar_products) const
{
#if HAVE_MKL || HAVE_CBLAS
XT::Common::Blas::dgemv(XT::Common::Blas::row_major(),
@@ -465,7 +482,7 @@ public:
#endif
}
- void apply_exponential(std::vector<RangeFieldType>& values) const
+ void apply_exponential(QuadratureWeightsType& values) const
{
assert(values.size() < std::numeric_limits<int>::max());
XT::Common::Mkl::exp(static_cast<int>(values.size()), values.data(), values.data());
@@ -571,10 +588,11 @@ public:
} // ii
} // void calculate_A_Binv(...)
- void apply_inverse_hessian(const DomainType& alpha, const DomainType& u, DomainType& Hinv_u) const
+ void
+ apply_inverse_hessian(const QuadratureWeightsType& density_evaluations, const DomainType& u, DomainType& Hinv_u) const
{
thread_local auto H = XT::Common::make_unique<MatrixType>();
- calculate_hessian(alpha, M_, *H);
+ calculate_hessian(density_evaluations, M_, *H);
XT::LA::cholesky(*H);
thread_local DomainType tmp_vec;
XT::LA::solve_lower_triangular(*H, tmp_vec, u);
@@ -600,11 +618,33 @@ public:
for (size_t ii = 0; ii < basis_dimRange; ++ii) {
auto* H_row = &(H[ii][0]);
const auto factor_ll_ii = basis_ll[ii] * factor_ll;
- if (!XT::Common::is_zero(factor_ll_ii)) {
- for (size_t kk = 0; kk <= ii; ++kk) {
- H_row[kk] += basis_ll[kk] * factor_ll_ii;
- } // kk
- }
+ // if (!XT::Common::is_zero(factor_ll_ii)) {
+ for (size_t kk = 0; kk <= ii; ++kk) {
+ H_row[kk] += basis_ll[kk] * factor_ll_ii;
+ } // kk
+ // }
+ } // ii
+ } // ll
+ } // void calculate_hessian(...)
+
+ void calculate_hessian(const QuadratureWeightsType& density_evaluations,
+ const BasisValuesMatrixType& M,
+ MatrixType& H) const
+ {
+ std::fill(H.begin(), H.end(), 0.);
+ const size_t num_quad_points = quad_weights_.size();
+ // matrix is symmetric, we only use lower triangular part
+ for (size_t ll = 0; ll < num_quad_points; ++ll) {
+ auto factor_ll = density_evaluations[ll] * quad_weights_[ll];
+ const auto* basis_ll = &(M.get_entry_ref(ll, 0.));
+ for (size_t ii = 0; ii < basis_dimRange; ++ii) {
+ auto* H_row = &(H[ii][0]);
+ const auto factor_ll_ii = basis_ll[ii] * factor_ll;
+ // if (!XT::Common::is_zero(factor_ll_ii)) {
+ for (size_t kk = 0; kk <= ii; ++kk) {
+ H_row[kk] += basis_ll[kk] * factor_ll_ii;
+ } // kk
+ // }
} // ii
} // ll
} // void calculate_hessian(...)
@@ -663,8 +703,8 @@ public:
} // void change_basis(...)
const MomentBasis& basis_functions_;
- std::vector<BasisDomainType> quad_points_;
- std::vector<RangeFieldType> quad_weights_;
+ QuadraturePointsType quad_points_;
+ QuadratureWeightsType quad_weights_;
BasisValuesMatrixType M_;
const RangeFieldType tau_;
const RangeFieldType epsilon_gamma_;
@@ -699,9 +739,10 @@ public:
using typename BaseType::BasisDomainType;
using typename BaseType::BasisValuesMatrixType;
using typename BaseType::DomainType;
+ using typename BaseType::FluxDomainType;
using typename BaseType::MatrixType;
using typename BaseType::MomentBasis;
- using typename BaseType::QuadraturePointsValuesType;
+ using typename BaseType::QuadratureWeightsType;
using typename BaseType::RangeFieldType;
using typename BaseType::VectorType;
@@ -732,43 +773,46 @@ public:
return ret;
}
- static RangeFieldType minmod(const RangeFieldType& first_slope, const RangeFieldType& second_slope)
+ DomainType get_u(const QuadratureWeightsType& density_evaluations) const
{
- 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;
+ DomainType ret(0.);
+ const size_t num_quad_points = quad_weights_.size();
+ for (size_t ll = 0; ll < num_quad_points; ++ll) {
+ const auto factor_ll = density_evaluations[ll] * quad_weights_[ll];
+ 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;
+ } // ll
return ret;
}
- static RangeFieldType maxmod(const RangeFieldType first_slope, const RangeFieldType second_slope)
+ void store_density_evaluations(QuadratureWeightsType& density_evaluations, const DomainType& alpha) const
{
- 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;
+ density_evaluations.resize(quad_points_.size());
+ this->calculate_scalar_products(alpha, M_, density_evaluations);
+ this->apply_exponential(density_evaluations);
}
- static RangeFieldType superbee(const RangeFieldType first_slope, const RangeFieldType second_slope)
+ void store_evaluations(QuadratureWeightsType& density_evaluations,
+ const std::function<RangeFieldType(const FluxDomainType&)>& boundary_density) const
{
- 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);
+ density_evaluations.resize(quad_points_.size());
+ for (size_t ll = 0; ll < quad_points_.size(); ++ll)
+ density_evaluations[ll] = boundary_density(quad_points_[ll]);
}
template <class IntersectionType>
- DomainType calculate_boundary_flux(const DomainType& alpha, const IntersectionType& intersection)
+ DomainType calculate_boundary_flux(const QuadratureWeightsType& density_evaluations,
+ 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;
+ RangeFieldType factor = density_evaluations[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;
@@ -777,38 +821,29 @@ public:
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
+ template <SlopeType slope_type, class FluxesMapType>
+ void calculate_reconstructed_fluxes(const FieldVector<const QuadratureWeightsType*, 3>& density_evaluations,
+ 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
+ thread_local XT::Common::FieldVector<QuadratureWeightsType, 2> reconstructed_values(
+ QuadratureWeightsType(quad_points_.size()));
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
+ if (slope_type == SlopeType::no_slope) {
+ for (size_t ll = 0; ll < quad_points_.size(); ++ll)
+ vals_left[ll] = vals_right[ll] = (*density_evaluations[1])[ll];
+ } else {
+ const auto slope_func =
+ (slope_type == SlopeType::minmod) ? XT::Common::minmod<RangeFieldType> : superbee<RangeFieldType>;
+ for (size_t ll = 0; ll < quad_points_.size(); ++ll) {
+ const auto slope = slope_func((*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;
@@ -826,9 +861,7 @@ public:
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(...)
+ } // void calculate_reconstructed_fluxes(...)
// returns (alpha, (actual_u, r)), where r is the regularization parameter and actual_u the regularized u
virtual std::unique_ptr<AlphaReturnType>
@@ -1221,10 +1254,9 @@ public:
using BasisValuesMatrixType = FieldVector<XT::LA::CommonDenseMatrix<RangeFieldType>, num_blocks>;
using QuadraturePointsType =
FieldVector<std::vector<BasisDomainType, boost::alignment::aligned_allocator<BasisDomainType, 64>>, num_blocks>;
- using QuadratureWeightsType =
- FieldVector<std::vector<RangeFieldType, boost::alignment::aligned_allocator<RangeFieldType, 64>>, num_blocks>;
- using TemporaryVectorType = std::vector<RangeFieldType, boost::alignment::aligned_allocator<RangeFieldType, 64>>;
- using TemporaryVectorsType = FieldVector<TemporaryVectorType, num_blocks>;
+ using BlockQuadratureWeightsType =
+ std::vector<RangeFieldType, boost::alignment::aligned_allocator<RangeFieldType, 64>>;
+ using QuadratureWeightsType = FieldVector<BlockQuadratureWeightsType, num_blocks>;
using AlphaReturnType = std::pair<BlockVectorType, std::pair<DomainType, RangeFieldType>>;
explicit EntropyBasedFluxImplementation(const MomentBasis& basis_functions,
@@ -1258,10 +1290,6 @@ public:
}
} // jj
for (size_t jj = 0; jj < num_blocks; ++jj) {
- while (quad_weights_[jj].size() % 8) { // align to 64 byte boundary
- quad_points_[jj].push_back(quad_points_[jj].back());
- quad_weights_[jj].push_back(0.);
- }
M_[jj] = XT::LA::CommonDenseMatrix<RangeFieldType>(quad_points_[jj].size(), block_size, 0., 0);
for (size_t ll = 0; ll < quad_points_[jj].size(); ++ll) {
const auto val = basis_functions_.evaluate(quad_points_[jj][ll], jj);
@@ -1338,7 +1366,7 @@ public:
const size_t dd) const
{
// calculate \sum_{i=1}^d < \omega_i m G_\alpha(u) > n_i
- thread_local FieldVector<TemporaryVectorsType, 2> work_vecs;
+ thread_local FieldVector<QuadratureWeightsType, 2> work_vecs;
for (size_t jj = 0; jj < num_blocks; ++jj) {
work_vecs[0][jj].resize(quad_points_[jj].size());
work_vecs[1][jj].resize(quad_points_[jj].size());
@@ -1378,6 +1406,108 @@ public:
return u_block;
}
+ DomainType get_u(const QuadratureWeightsType& density_evaluations) const
+ {
+ DomainType ret;
+ for (size_t jj = 0; jj < num_blocks; ++jj) {
+ const auto offset = block_size * jj;
+ for (size_t ll = 0; ll < quad_weights_[jj].size(); ++ll) {
+ const auto factor = density_evaluations[jj][ll] * quad_weights_[jj][ll];
+ const auto* basis_ll = &(M_[jj].get_entry_ref(ll, 0.));
+ for (size_t ii = 0; ii < block_size; ++ii)
+ ret[offset + ii] += basis_ll[ii] * factor;
+ } // ll
+ } // jj (intervals)
+ return ret;
+ }
+
+ void store_density_evaluations(QuadratureWeightsType& density_evaluations, const DomainType& alpha) const
+ {
+ this->calculate_scalar_products(alpha, M_, density_evaluations);
+ this->apply_exponential(density_evaluations);
+ }
+
+ void store_evaluations(QuadratureWeightsType& density_evaluations,
+ const std::function<RangeFieldType(const FluxDomainType&)>& boundary_density) const
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj) {
+ density_evaluations[jj].resize(quad_points_[jj].size());
+ for (size_t ll = 0; ll < quad_points_[jj].size(); ++ll)
+ density_evaluations[jj][ll] = boundary_density(quad_points_[jj][ll]);
+ } // jj
+ }
+
+ template <class IntersectionType>
+ DomainType calculate_boundary_flux(const QuadratureWeightsType& density_evaluations,
+ 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_blocks; ++jj) {
+ const auto offset = block_size * jj;
+ for (size_t ll = 0; ll < quad_points_[jj].size(); ++ll) {
+ const auto position = quad_points_[jj][ll][dd];
+ if (position * n_ij[dd] < 0.) {
+ RangeFieldType factor = density_evaluations[jj][ll] * quad_weights_[jj][ll] * position;
+ for (size_t ii = 0; ii < block_size; ++ii)
+ ret[offset + ii] += M_[jj].get_entry(ll, ii) * factor;
+ }
+ } // ll
+ } // jj
+ return ret;
+ }
+
+ template <SlopeType slope_type, class FluxesMapType>
+ void calculate_reconstructed_fluxes(const FieldVector<const QuadratureWeightsType*, 3>& density_evaluations,
+ FluxesMapType& flux_values,
+ const size_t dd) const
+ {
+ // 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.);
+
+ // evaluate exp(alpha^T b(v_i)) at all quadratures points v_i for all three alphas
+ thread_local XT::Common::FieldVector<BlockQuadratureWeightsType, 2> reconstructed_values(
+ BlockQuadratureWeightsType(quad_points_[0].size()));
+
+ // get left and right reconstructed values for each quadrature point v_i
+ auto& vals_left = reconstructed_values[0];
+ auto& vals_right = reconstructed_values[1];
+
+ const auto slope_func =
+ (slope_type == SlopeType::minmod) ? XT::Common::minmod<RangeFieldType> : superbee<RangeFieldType>;
+ for (size_t jj = 0; jj < num_blocks; ++jj) {
+ // reconstruct densities
+ if (slope_type == SlopeType::no_slope) {
+ for (size_t ll = 0; ll < quad_points_[jj].size(); ++ll)
+ vals_left[ll] = vals_right[ll] = (*density_evaluations[1])[jj][ll];
+ } else {
+ for (size_t ll = 0; ll < quad_points_[jj].size(); ++ll) {
+ const auto slope = slope_func((*density_evaluations[1])[jj][ll] - (*density_evaluations[0])[jj][ll],
+ (*density_evaluations[2])[jj][ll] - (*density_evaluations[1])[jj][ll]);
+ vals_left[ll] = (*density_evaluations[1])[jj][ll] - 0.5 * slope;
+ vals_right[ll] = (*density_evaluations[1])[jj][ll] + 0.5 * slope;
+ } // ll
+ }
+ // calculate fluxes
+ const auto offset = block_size * jj;
+ for (size_t ll = 0; ll < quad_points_[jj].size(); ++ll) {
+ const auto position = quad_points_[jj][ll][dd];
+ 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;
+ for (size_t ii = 0; ii < block_size; ++ii)
+ val[offset + ii] += M_[jj].get_entry(ll, ii) * factor;
+ } // ll
+ } // jj
+ } // void calculate_reconstructed_fluxes(...)
+
std::unique_ptr<AlphaReturnType>
get_alpha(const DomainType& u, const VectorType& alpha_in, const bool regularize) const
{
@@ -1518,9 +1648,9 @@ public:
}
// temporary vectors to store inner products and exponentials
- TemporaryVectorsType& working_storage() const
+ QuadratureWeightsType& working_storage() const
{
- thread_local TemporaryVectorsType work_vecs;
+ thread_local QuadratureWeightsType work_vecs;
for (size_t jj = 0; jj < num_blocks; ++jj)
work_vecs[jj].resize(quad_points_[jj].size());
return work_vecs;
@@ -1542,7 +1672,7 @@ public:
void calculate_scalar_products_block(const size_t jj,
const LocalVectorType& beta_in,
const XT::LA::CommonDenseMatrix<RangeFieldType>& M,
- TemporaryVectorType& scalar_products) const
+ BlockQuadratureWeightsType& scalar_products) const
{
const size_t num_quad_points = quad_points_[jj].size();
for (size_t ll = 0; ll < num_quad_points; ++ll) {
@@ -1553,19 +1683,21 @@ public:
void calculate_scalar_products(const BlockVectorType& beta_in,
const BasisValuesMatrixType& M,
- TemporaryVectorsType& scalar_products) const
+ QuadratureWeightsType& scalar_products) const
{
- for (size_t jj = 0; jj < num_blocks; ++jj)
+ for (size_t jj = 0; jj < num_blocks; ++jj) {
+ scalar_products[jj].resize(quad_weights_[jj].size());
calculate_scalar_products_block(jj, beta_in.block(jj), M[jj], scalar_products[jj]);
+ }
}
- void apply_exponential(TemporaryVectorType& values) const
+ void apply_exponential(BlockQuadratureWeightsType& values) const
{
assert(values.size() < std::numeric_limits<int>::max());
XT::Common::Mkl::exp(static_cast<int>(values.size()), values.data(), values.data());
}
- void apply_exponential(TemporaryVectorsType& values) const
+ void apply_exponential(QuadratureWeightsType& values) const
{
for (size_t jj = 0; jj < num_blocks; ++jj)
apply_exponential(values[jj]);
@@ -1772,10 +1904,11 @@ public:
} // jj
} // void calculate_A_Binv(...)
- void apply_inverse_hessian(const DomainType& alpha, const DomainType& u, DomainType& Hinv_u) const
+ void
+ apply_inverse_hessian(const QuadratureWeightsType& density_evaluations, const DomainType& u, DomainType& Hinv_u) const
{
thread_local auto H = XT::Common::make_unique<BlockMatrixType>();
- calculate_hessian(alpha, M_, *H);
+ calculate_hessian(density_evaluations, M_, *H);
for (size_t jj = 0; jj < num_blocks; ++jj)
XT::LA::cholesky(H->block(jj));
FieldVector<RangeFieldType, block_size> tmp_vec;
@@ -1792,17 +1925,16 @@ public:
} // jj
}
- void calculate_hessian(const BlockVectorType& alpha, const BasisValuesMatrixType& M, BlockMatrixType& H) const
+ void calculate_hessian(const QuadratureWeightsType& density_evaluations,
+ const BasisValuesMatrixType& M,
+ BlockMatrixType& H) const
{
- auto& work_vec = working_storage();
- calculate_scalar_products(alpha, M, work_vec);
- apply_exponential(work_vec);
// matrix is symmetric, we only use lower triangular part
for (size_t jj = 0; jj < num_blocks; ++jj) {
std::fill(H.block(jj).begin(), H.block(jj).end(), 0.);
const size_t num_quad_points = quad_weights_[jj].size();
for (size_t ll = 0; ll < num_quad_points; ++ll) {
- auto factor_ll = work_vec[jj][ll] * quad_weights_[jj][ll];
+ auto factor_ll = density_evaluations[jj][ll] * quad_weights_[jj][ll];
const auto* basis_ll = &(M[jj].get_entry_ref(ll, 0.));
for (size_t ii = 0; ii < block_size; ++ii) {
auto* H_row = &(H.block(jj)[ii][0]);
@@ -1814,6 +1946,14 @@ public:
} // jj
} // void calculate_hessian(...)
+ void calculate_hessian(const BlockVectorType& alpha, const BasisValuesMatrixType& M, BlockMatrixType& H) const
+ {
+ auto& work_vec = working_storage();
+ calculate_scalar_products(alpha, M, work_vec);
+ apply_exponential(work_vec);
+ calculate_hessian(work_vec, M, H);
+ } // void calculate_hessian(...)
+
// J = df/dalpha is the derivative of the flux with respect to alpha.
// As F = (f_1, f_2, f_3) is matrix-valued
// (div f = \sum_{i=1}^d \partial_{x_i} f_i = \sum_{i=1}^d \partial_{x_i} < v_i m \hat{psi}(alpha) > is
@@ -2156,6 +2296,135 @@ public:
return u;
}
+ DomainType get_u(const QuadratureWeightsType& density_evaluations) const
+ {
+ DomainType u(0.);
+ LocalVectorType local_u;
+ const auto& triangulation = basis_functions_.triangulation();
+ const auto& faces = triangulation.faces();
+ for (size_t jj = 0; jj < faces.size(); ++jj) {
+ const auto& face = faces[jj];
+ const auto& vertices = face->vertices();
+ local_u *= 0.;
+ for (size_t ll = 0; ll < quad_weights_[jj].size(); ++ll) {
+ const auto& basis_ll = M_[jj][ll];
+ for (size_t ii = 0; ii < 3; ++ii)
+ local_u[ii] += basis_ll[ii] * density_evaluations[jj][ll];
+ } // ll (quad points)
+ for (size_t ii = 0; ii < 3; ++ii)
+ u[vertices[ii]->index()] = local_u[ii];
+ } // jj (faces)
+ return u;
+ }
+
+ void store_density_evaluations(QuadratureWeightsType& density_evaluations, const DomainType& alpha) const
+ {
+ const auto& triangulation = basis_functions_.triangulation();
+ const auto& faces = triangulation.faces();
+ density_evaluations.resize(faces.size());
+ XT::Common::FieldVector<RangeFieldType, 3> local_alpha;
+ for (size_t jj = 0; jj < faces.size(); ++jj) {
+ density_evaluations[jj].resize(quad_weights_[jj].size());
+ const auto& face = faces[jj];
+ const auto& vertices = face->vertices();
+ for (size_t ii = 0; ii < 3; ++ii)
+ local_alpha[ii] = alpha[vertices[ii]->index()];
+ for (size_t ll = 0; ll < quad_weights_[jj].size(); ++ll)
+ density_evaluations[jj][ll] = std::exp(local_alpha * M_[jj][ll]);
+ } // jj (faces)
+ }
+
+ void store_evaluations(QuadratureWeightsType& density_evaluations,
+ const std::function<RangeFieldType(const FluxDomainType&)>& boundary_density) const
+ {
+ const auto& triangulation = basis_functions_.triangulation();
+ const auto& faces = triangulation.faces();
+ density_evaluations.resize(faces.size());
+ for (size_t jj = 0; jj < faces.size(); ++jj) {
+ density_evaluations[jj].resize(quad_points_[jj].size());
+ for (size_t ll = 0; ll < quad_points_[jj].size(); ++ll)
+ density_evaluations[jj][ll] = boundary_density(quad_points_[jj][ll]);
+ } // jj
+ }
+
+ template <class IntersectionType>
+ DomainType calculate_boundary_flux(const QuadratureWeightsType& density_evaluations,
+ 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();
+ const auto& triangulation = basis_functions_.triangulation();
+ const auto& faces = triangulation.faces();
+ DomainType ret(0.);
+ LocalVectorType local_ret;
+ for (size_t jj = 0; jj < faces.size(); ++jj) {
+ local_ret *= 0.;
+ const auto& face = faces[jj];
+ const auto& vertices = face->vertices();
+ for (size_t ll = 0; ll < quad_weights_[jj].size(); ++ll) {
+ const auto position = quad_points_[jj][ll][dd];
+ if (position * n_ij[dd] < 0.) {
+ const auto& basis_ll = M_[jj][ll];
+ RangeFieldType factor = density_evaluations[jj][ll] * quad_weights_[jj][ll] * position;
+ for (size_t ii = 0; ii < 3; ++ii)
+ local_ret[ii] += basis_ll[ii] * factor;
+ }
+ } // ll (quad points)
+ for (size_t ii = 0; ii < 3; ++ii)
+ ret[vertices[ii]->index()] += local_ret[ii];
+ } // jj
+ return ret;
+ }
+
+ template <SlopeType slope_type, class FluxesMapType>
+ void calculate_reconstructed_fluxes(const FieldVector<const QuadratureWeightsType*, 3>& density_evaluations,
+ FluxesMapType& flux_values,
+ const size_t dd) const
+ {
+ // 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.);
+ thread_local XT::Common::FieldVector<std::vector<RangeFieldType>, 2> reconstructed_values(
+ std::vector<RangeFieldType>(quad_points_[0].size()));
+ const auto& triangulation = basis_functions_.triangulation();
+ const auto& faces = triangulation.faces();
+ const auto slope_func =
+ (slope_type == SlopeType::minmod) ? XT::Common::minmod<RangeFieldType> : superbee<RangeFieldType>;
+ auto& vals_left = reconstructed_values[0];
+ auto& vals_right = reconstructed_values[1];
+ for (size_t jj = 0; jj < faces.size(); ++jj) {
+ const auto& face = faces[jj];
+ const auto& vertices = face->vertices();
+ // reconstruct densities
+ if (slope_type == SlopeType::no_slope) {
+ for (size_t ll = 0; ll < quad_weights_[jj].size(); ++ll)
+ vals_left[ll] = vals_right[ll] = (*density_evaluations[1])[jj][ll];
+ } else {
+ for (size_t ll = 0; ll < quad_weights_[jj].size(); ++ll) {
+ const auto slope = slope_func((*density_evaluations[1])[jj][ll] - (*density_evaluations[0])[jj][ll],
+ (*density_evaluations[2])[jj][ll] - (*density_evaluations[1])[jj][ll]);
+ vals_left[ll] = (*density_evaluations[1])[jj][ll] - 0.5 * slope;
+ vals_right[ll] = (*density_evaluations[1])[jj][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][dd];
+ 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 ii = 0; ii < 3; ++ii)
+ val[vertices[ii]->index()] += basis_ll[ii] * factor;
+ } // ll (quad points)
+ } // jj
+ } // void calculate_reconstructed_fluxes(...)
+
std::unique_ptr<AlphaReturnType>
get_alpha(const DomainType& u, const VectorType& alpha_in, const bool regularize) const
{
@@ -2301,7 +2570,7 @@ public:
}
// temporary vectors to store inner products and exponentials
- std::vector<std::vector<RangeFieldType>>& get_work_vecs() const
+ std::vector<std::vector<RangeFieldType>>& working_storage() const
{
thread_local std::vector<std::vector<RangeFieldType>> work_vecs;
const auto& triangulation = basis_functions_.triangulation();
@@ -2343,15 +2612,14 @@ public:
} // ii
} // void calculate_J_Hinv(...)
- void apply_inverse_hessian(const DomainType& alpha, const DomainType& u, DomainType& Hinv_u) const
+ void
+ apply_inverse_hessian(const QuadratureWeightsType& density_evaluations, 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);
+ calculate_hessian(density_evaluations, M_, H);
# if HAVE_EIGEN
+ thread_local VectorType u_vec(basis_dimRange, 0., 0), Hinv_u_vec(basis_dimRange, 0., 0);
+ std::copy(u.begin(), u.end(), u_vec.begin());
typedef ::Eigen::SparseMatrix<RangeFieldType, ::Eigen::ColMajor> ColMajorBackendType;
ColMajorBackendType colmajor_copy(H.backend());
typedef ::Eigen::SimplicialLDLT<ColMajorBackendType> SolverType;
@@ -2390,7 +2658,7 @@ public:
LocalVectorType local_alpha, local_u;
const auto& triangulation = basis_functions_.triangulation();
const auto& faces = triangulation.faces();
- auto& work_vecs = get_work_vecs();
+ auto& work_vecs = working_storage();
for (size_t jj = 0; jj < faces.size(); ++jj) {
const auto& face = faces[jj];
const auto& vertices = face->vertices();
@@ -2414,6 +2682,30 @@ public:
g_k -= v;
}
+ void calculate_hessian(const QuadratureWeightsType& density_evaluations,
+ const BasisValuesMatrixType& M,
+ SparseMatrixType& H) const
+ {
+ H *= 0.;
+ LocalMatrixType H_local(0.);
+ const auto& triangulation = basis_functions_.triangulation();
+ const auto& faces = triangulation.faces();
+ for (size_t jj = 0; jj < faces.size(); ++jj) {
+ H_local *= 0.;
+ const auto& face = faces[jj];
+ const auto& vertices = face->vertices();
+ for (size_t ll = 0; ll < quad_weights_[jj].size(); ++ll) {
+ const auto& basis_ll = M[jj][ll];
+ for (size_t ii = 0; ii < 3; ++ii)
+ for (size_t kk = 0; kk < 3; ++kk)
+ H_local[ii][kk] += basis_ll[ii] * basis_ll[kk] * density_evaluations[jj][ll] * quad_weights_[jj][ll];
+ } // ll (quad points)
+ for (size_t ii = 0; ii < 3; ++ii)
+ for (size_t kk = 0; kk < 3; ++kk)
+ H.add_to_entry(vertices[ii]->index(), vertices[kk]->index(), H_local[ii][kk]);
+ } // jj (faces)
+ } // void calculate_hessian(...)
+
void calculate_hessian(const VectorType& alpha,
const BasisValuesMatrixType& M,
SparseMatrixType& H,
@@ -2424,7 +2716,7 @@ public:
LocalMatrixType H_local(0.);
const auto& triangulation = basis_functions_.triangulation();
const auto& faces = triangulation.faces();
- auto& work_vecs = get_work_vecs();
+ auto& work_vecs = working_storage();
for (size_t jj = 0; jj < faces.size(); ++jj) {
H_local *= 0.;
const auto& face = faces[jj];
@@ -2457,7 +2749,7 @@ public:
assert(dd < dimFlux);
J_dd *= 0.;
LocalMatrixType J_local(0.);
- auto& work_vecs = get_work_vecs();
+ auto& work_vecs = working_storage();
const auto& triangulation = basis_functions_.triangulation();
const auto& faces = triangulation.faces();
for (size_t jj = 0; jj < faces.size(); ++jj) {
@@ -3270,7 +3562,7 @@ public:
using LocalVectorType = XT::Common::FieldVector<RangeFieldType, block_size>;
using BasisValuesMatrixType = FieldVector<std::vector<LocalVectorType>, num_intervals>;
using QuadraturePointsType = FieldVector<std::vector<RangeFieldType>, num_intervals>;
- using QuadratureWeightsType = FieldVector<std::vector<RangeFieldType>, num_intervals>;
+ using QuadratureWeightsType = QuadraturePointsType;
explicit EntropyBasedFluxImplementation(const MomentBasis& basis_functions,
const RangeFieldType tau,
@@ -3414,31 +3706,42 @@ public:
return ret;
}
- static RangeFieldType minmod(const RangeFieldType& first_slope, const RangeFieldType& second_slope)
+ DomainType get_u(const QuadratureWeightsType& density_evaluations) const
{
- 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;
+ DomainType ret;
+ for (size_t jj = 0; jj < num_intervals; ++jj) {
+ for (size_t ll = 0; ll < quad_weights_[jj].size(); ++ll) {
+ const auto& basis_ll = M_[jj][ll];
+ auto factor_ll = density_evaluations[jj][ll] * quad_weights_[jj][ll];
+ for (size_t ii = 0; ii < 2; ++ii)
+ ret[jj + ii] += basis_ll[ii] * factor_ll;
+ } // ll (quad points)
+ } // jj (intervals)
return ret;
}
- static RangeFieldType maxmod(const RangeFieldType first_slope, const RangeFieldType second_slope)
+ void store_density_evaluations(QuadratureWeightsType& density_evaluations, const DomainType& alpha) const
{
- 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;
+ for (size_t jj = 0; jj < num_intervals; ++jj) {
+ density_evaluations[jj].resize(quad_points_[jj].size());
+ for (size_t ll = 0; ll < quad_points_[jj].size(); ++ll)
+ density_evaluations[jj][ll] = std::exp(alpha[jj] * M_[jj][ll][0] + alpha[jj + 1] * M_[jj][ll][1]);
+ } // jj
}
- static RangeFieldType superbee(const RangeFieldType first_slope, const RangeFieldType second_slope)
+ void store_evaluations(QuadratureWeightsType& density_evaluations,
+ const std::function<RangeFieldType(const FluxDomainType&)>& boundary_density) const
{
- 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);
+ for (size_t jj = 0; jj < num_intervals; ++jj) {
+ density_evaluations[jj].resize(quad_points_[jj].size());
+ for (size_t ll = 0; ll < quad_points_[jj].size(); ++ll)
+ density_evaluations[jj][ll] = boundary_density(quad_points_[jj][ll]);
+ } // jj
}
template <class IntersectionType>
- DomainType calculate_boundary_flux(const DomainType& alpha, const IntersectionType& intersection)
+ DomainType calculate_boundary_flux(const QuadratureWeightsType& density_evaluations,
+ const IntersectionType& intersection)
{
// evaluate exp(alpha^T b(v_i)) at all quadratures points v_i
const auto dd = intersection.indexInInside() / 2;
@@ -3448,8 +3751,7 @@ public:
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;
+ RangeFieldType factor = density_evaluations[jj][ll] * 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;
@@ -3459,18 +3761,14 @@ public:
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
+ template <SlopeType slope_type, class FluxesMapType>
+ void calculate_reconstructed_fluxes(const FieldVector<const QuadratureWeightsType*, 3>& density_evaluations,
+ 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);
@@ -3480,29 +3778,23 @@ public:
auto& right_flux_value = flux_values[coord];
right_flux_value = left_flux_value = DomainType(0.);
+ auto& vals_left = reconstructed_values[0];
+ auto& vals_right = reconstructed_values[1];
+ const auto slope_func =
+ (slope_type == SlopeType::minmod) ? XT::Common::minmod<RangeFieldType> : superbee<RangeFieldType>;
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)
-
+ if (slope_type == SlopeType::no_slope) {
+ for (size_t ll = 0; ll < quad_weights_[jj].size(); ++ll)
+ vals_left[ll] = vals_right[ll] = (*density_evaluations[1])[jj][ll];
+ } else {
+ for (size_t ll = 0; ll < quad_weights_[jj].size(); ++ll) {
+ const auto slope = slope_func((*density_evaluations[1])[jj][ll] - (*density_evaluations[0])[jj][ll],
+ (*density_evaluations[2])[jj][ll] - (*density_evaluations[1])[jj][ll]);
+ vals_left[ll] = (*density_evaluations[1])[jj][ll] - 0.5 * slope;
+ vals_right[ll] = (*density_evaluations[1])[jj][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];
@@ -3514,7 +3806,7 @@ public:
val[jj + mm] += basis_ll[mm] * factor;
} // ll (quad points)
} // jj (intervals)
- } // void calculate_minmod_density_reconstruction(...)
+ } // void calculate_reconstructed_fluxes(...)
// returns (alpha, (actual_u, r)), where r is the regularization parameter and actual_u the regularized u
std::unique_ptr<AlphaReturnType>
@@ -3634,15 +3926,14 @@ public:
return true;
}
- FieldVector<std::vector<RangeFieldType>, num_intervals>& working_storage() const
+ QuadratureWeightsType& working_storage() const
{
- thread_local FieldVector<std::vector<RangeFieldType>, num_intervals> work_vec;
+ thread_local QuadratureWeightsType work_vec;
for (size_t jj = 0; jj < num_intervals; ++jj)
work_vec[jj].resize(quad_points_[jj].size());
return work_vec;
}
-
RangeFieldType calculate_f(const VectorType& alpha, const VectorType& v) const
{
RangeFieldType ret(0.);
@@ -3679,6 +3970,24 @@ public:
g_k -= v;
}
+ void calculate_hessian(const QuadratureWeightsType& density_evalutations,
+ const BasisValuesMatrixType& M,
+ VectorType& H_diag,
+ FieldVector<RangeFieldType, dimRange - 1>& H_subdiag) const
+ {
+ std::fill(H_diag.begin(), H_diag.end(), 0.);
+ std::fill(H_subdiag.begin(), H_subdiag.end(), 0.);
+ for (size_t jj = 0; jj < num_intervals; ++jj) {
+ for (size_t ll = 0; ll < quad_weights_[jj].size(); ++ll) {
+ const auto& basis_ll = M[jj][ll];
+ const auto factor = density_evalutations[jj][ll] * quad_weights_[jj][ll];
+ for (size_t ii = 0; ii < 2; ++ii)
+ H_diag[jj + ii] += std::pow(basis_ll[ii], 2) * factor;
+ H_subdiag[jj] += basis_ll[0] * basis_ll[1] * factor;
+ } // ll (quad points)
+ } // jj (intervals)
+ } // void calculate_hessian(...)
+
void calculate_hessian(const DomainType& alpha,
const BasisValuesMatrixType& M,
VectorType& H_diag,
@@ -3725,12 +4034,12 @@ public:
} // jj (intervals)
} // void calculate_J(...)
- void apply_inverse_hessian(const DomainType& alpha, const DomainType& u, DomainType& Hinv_u) const
+ void
+ apply_inverse_hessian(const QuadratureWeightsType& density_evaluations, 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;
+ calculate_hessian(density_evaluations, M_, H_diag, H_subdiag);
// 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
diff --git a/dune/gdt/momentmodels/entropyflux_kineticcoords.hh b/dune/gdt/momentmodels/entropyflux_kineticcoords.hh
index bb65eab97..889b05e6d 100644
--- a/dune/gdt/momentmodels/entropyflux_kineticcoords.hh
+++ b/dune/gdt/momentmodels/entropyflux_kineticcoords.hh
@@ -26,7 +26,7 @@ namespace Dune {
namespace GDT {
-template <class GridViewImp, class MomentBasisImp, bool fluxes_precomputed>
+template <class GridViewImp, class MomentBasisImp, SlopeType slope>
class EntropyBasedFluxEntropyCoordsFunction
: public XT::Functions::FluxFunctionInterface<XT::Grid::extract_entity_t<GridViewImp>,
MomentBasisImp::dimRange,
@@ -56,6 +56,9 @@ public:
using AlphaReturnType = typename ImplementationType::AlphaReturnType;
using VectorType = typename ImplementationType::VectorType;
using I = XT::Grid::extract_intersection_t<GridViewType>;
+ using QuadratureWeightsType = typename ImplementationType::QuadratureWeightsType;
+ using BoundaryQuadratureWeightsType =
+ std::vector<XT::Common::FieldVector<XT::Common::FieldVector<QuadratureWeightsType, 2>, dimFlux>>;
explicit EntropyBasedFluxEntropyCoordsFunction(
const MomentBasis& basis_functions,
@@ -75,7 +78,6 @@ public:
: implementation_(other.implementation_)
{}
-
static const constexpr bool available = true;
class Localfunction : public LocalFunctionType
@@ -132,17 +134,12 @@ public:
StateType evaluate_kinetic_flux(const E& /*inside_entity*/,
const E& /*outside_entity*/,
- // const StateType& alpha_i,
- // const StateType& alpha_j,
- const StateType& flux_or_alpha_i,
- const StateType& flux_or_alpha_j,
+ const StateType& flux_i,
+ const StateType& flux_j,
const DomainType& n_ij,
const size_t dd) const
{
- if (fluxes_precomputed)
- return (flux_or_alpha_i + flux_or_alpha_j) * n_ij[dd];
- else
- return implementation_->evaluate_kinetic_flux_with_alphas(flux_or_alpha_i, flux_or_alpha_j, n_ij, dd);
+ return (flux_i + flux_j) * n_ij[dd];
} // StateType evaluate_kinetic_flux(...)
const MomentBasis& basis_functions() const
@@ -155,6 +152,11 @@ public:
return implementation_->get_u(alpha);
}
+ StateType get_u(const size_t entity_index) const
+ {
+ return implementation_->get_u(density_evaluations_[entity_index]);
+ }
+
StateType get_alpha(const StateType& u) const
{
const auto alpha = implementation_->get_alpha(u)->first;
@@ -164,25 +166,71 @@ public:
}
template <class FluxesMapType>
- void calculate_minmod_density_reconstruction(const FieldVector<StateType, 3>& alphas,
- FluxesMapType& precomputed_fluxes,
- const size_t dd) const
+ void calculate_reconstructed_fluxes(const FieldVector<size_t, 3>& entity_indices,
+ const FieldVector<bool, 3>& boundary_direction,
+ FluxesMapType& precomputed_fluxes,
+ const size_t dd) const
+ {
+ FieldVector<const QuadratureWeightsType*, 3> densities_stencil;
+ for (size_t ii = 0; ii < 3; ++ii)
+ if (entity_indices[ii] == size_t(-1))
+ densities_stencil[ii] = &boundary_density_evaluations_[entity_indices[1]][dd][boundary_direction[ii]];
+ else
+ densities_stencil[ii] = &density_evaluations_[entity_indices[ii]];
+ implementation_->template calculate_reconstructed_fluxes<slope, FluxesMapType>(
+ densities_stencil, precomputed_fluxes, dd);
+ }
+
+ StateType calculate_boundary_flux(const size_t entity_index, const I& intersection)
+ {
+ const size_t dd = intersection.indexInInside() / 2;
+ const size_t dir = intersection.indexInInside() % 2;
+ return implementation_->calculate_boundary_flux(boundary_density_evaluations_[entity_index][dd][dir], intersection);
+ }
+
+ void apply_inverse_hessian(const size_t entity_index, const StateType& u, StateType& Hinv_u) const
{
- implementation_->calculate_minmod_density_reconstruction(alphas, precomputed_fluxes, dd);
+ implementation_->apply_inverse_hessian(density_evaluations_[entity_index], u, Hinv_u);
}
- StateType calculate_boundary_flux(const StateType& alpha, const I& intersection)
+ void store_density_evaluations(const size_t entity_index, const StateType& alpha)
{
- return implementation_->calculate_boundary_flux(alpha, intersection);
+ implementation_->store_density_evaluations(density_evaluations_[entity_index], alpha);
}
- void apply_inverse_hessian(const StateType& alpha, const StateType& u, StateType& Hinv_u) const
+ void store_boundary_evaluations(const std::function<RangeFieldType(const DomainType&)>& boundary_density,
+ const size_t entity_index,
+ const size_t intersection_index)
{
- implementation_->apply_inverse_hessian(alpha, u, Hinv_u);
+ implementation_->store_evaluations(
+ boundary_density_evaluations_[entity_index][intersection_index / 2][intersection_index % 2], boundary_density);
}
+ std::vector<QuadratureWeightsType>& density_evaluations()
+ {
+ return density_evaluations_;
+ }
+
+ const std::vector<QuadratureWeightsType>& density_evaluations() const
+ {
+ return density_evaluations_;
+ }
+
+ BoundaryQuadratureWeightsType& boundary_density_evaluations()
+ {
+ return boundary_density_evaluations_;
+ }
+
+ const BoundaryQuadratureWeightsType& boundary_density_evaluations() const
+ {
+ return boundary_density_evaluations_;
+ }
+
+
private:
std::shared_ptr<ImplementationType> implementation_;
+ std::vector<QuadratureWeightsType> density_evaluations_;
+ BoundaryQuadratureWeightsType boundary_density_evaluations_;
};
diff --git a/dune/gdt/momentmodels/hessianinverter.hh b/dune/gdt/momentmodels/hessianinverter.hh
index 3f0320516..b3c525939 100644
--- a/dune/gdt/momentmodels/hessianinverter.hh
+++ b/dune/gdt/momentmodels/hessianinverter.hh
@@ -24,16 +24,15 @@ namespace Dune {
namespace GDT {
-template <class SpaceType, class VectorType, class MomentBasis, bool reconstruction>
+template <class SpaceType, class VectorType, class MomentBasis, SlopeType slope>
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, reconstruction>;
+ using EntropyFluxType = EntropyBasedFluxEntropyCoordsFunction<GridViewType, MomentBasis, slope>;
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>;
@@ -45,9 +44,7 @@ public:
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 = "")
+ const XT::Common::Parameter& param)
: space_(space)
, alpha_(space_, alpha_dofs, "alpha")
, u_update_(space_, u_update_dofs, "u_update")
@@ -56,11 +53,7 @@ public:
, 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)
@@ -73,11 +66,7 @@ public:
, 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
@@ -87,15 +76,12 @@ public:
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) {
+ XT::Common::FieldVector<RangeFieldType, dimRange> u, Hinv_u;
+ 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);
+ analytical_flux_.apply_inverse_hessian(space_.grid_view().indexSet().index(entity), 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: "
@@ -116,16 +102,12 @@ private:
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,
- bool reconstruction,
+ SlopeType slope,
class MatrixType = typename XT::LA::Container<typename MomentBasisImp::RangeFieldType>::MatrixType>
class EntropicHessianInverter
: public OperatorInterface<MatrixType, typename SpaceImp::GridViewType, MomentBasisImp::dimRange, 1>
@@ -138,19 +120,12 @@ public:
using SpaceType = SpaceImp;
using SourceSpaceType = SpaceImp;
using RangeSpaceType = SpaceImp;
- using EntropyFluxType =
- EntropyBasedFluxEntropyCoordsFunction<typename SpaceType::GridViewType, MomentBasis, reconstruction>;
+ using EntropyFluxType = EntropyBasedFluxEntropyCoordsFunction<typename SpaceType::GridViewType, MomentBasis, slope>;
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 = "")
+ EntropicHessianInverter(const EntropyFluxType& analytical_flux, const SpaceType& space)
: analytical_flux_(analytical_flux)
, space_(space)
- , min_acceptable_density_(min_acceptable_density)
- , filename_(filename)
{}
virtual bool linear() const override final
@@ -180,8 +155,8 @@ public:
VectorType& alpha_update,
const XT::Common::Parameter& param) const
{
- LocalEntropicHessianInverter<SpaceType, VectorType, MomentBasis, reconstruction> local_hessian_inverter(
- space_, alpha, u_update, alpha_update, analytical_flux_, min_acceptable_density_, param, filename_);
+ LocalEntropicHessianInverter<SpaceType, VectorType, MomentBasis, slope> local_hessian_inverter(
+ space_, alpha, u_update, alpha_update, analytical_flux_, param);
auto walker = XT::Grid::Walker<typename SpaceType::GridViewType>(space_.grid_view());
walker.append(local_hessian_inverter);
walker.walk(true);
@@ -190,8 +165,6 @@ public:
private:
const EntropyFluxType& analytical_flux_;
const SpaceType& space_;
- const RangeFieldType min_acceptable_density_;
- const std::string filename_;
}; // class EntropicHessianInverter<...>
diff --git a/dune/gdt/operators/advection-fv-entropybased.hh b/dune/gdt/operators/advection-fv-entropybased.hh
index e9a5b6fe1..06f37a1bb 100644
--- a/dune/gdt/operators/advection-fv-entropybased.hh
+++ b/dune/gdt/operators/advection-fv-entropybased.hh
@@ -64,7 +64,7 @@ public:
const InverseHessianOperatorType& inverse_hessian_operator_;
}; // class EntropicCoordinatesOperator<...>
-template <class AdvectionOperatorImp, class RhsOperatorImp, class InverseHessianOperatorImp>
+template <class DensityOperatorImp, class AdvectionOperatorImp, class RhsOperatorImp, class InverseHessianOperatorImp>
class EntropicCoordinatesCombinedOperator
: public OperatorInterface<typename AdvectionOperatorImp::MatrixType,
typename AdvectionOperatorImp::SGV,
@@ -79,14 +79,17 @@ public:
using typename BaseType::SourceSpaceType;
using typename BaseType::VectorType;
+ using DensityOperatorType = DensityOperatorImp;
using AdvectionOperatorType = AdvectionOperatorImp;
using RhsOperatorType = RhsOperatorImp;
using InverseHessianOperatorType = InverseHessianOperatorImp;
- EntropicCoordinatesCombinedOperator(const AdvectionOperatorType& advection_op,
+ EntropicCoordinatesCombinedOperator(const DensityOperatorType& density_op,
+ const AdvectionOperatorType& advection_op,
const RhsOperatorType& rhs_op,
const InverseHessianOperatorType& inverse_hessian_operator)
- : advection_op_(advection_op)
+ : density_op_(density_op)
+ , advection_op_(advection_op)
, rhs_op_(rhs_op)
, inverse_hessian_operator_(inverse_hessian_operator)
{}
@@ -108,6 +111,7 @@ public:
void apply(const VectorType& source, VectorType& range, const XT::Common::Parameter& param) const override final
{
+ density_op_.apply(source, range, param);
VectorType u_update = range;
std::fill(u_update.begin(), u_update.end(), 0.);
advection_op_.apply(source, u_update, param);
@@ -116,6 +120,7 @@ public:
inverse_hessian_operator_.apply_inverse_hessian(source, u_update, range, param);
}
+ const DensityOperatorType& density_op_;
const AdvectionOperatorType& advection_op_;
const RhsOperatorType& rhs_op_;
const InverseHessianOperatorType& inverse_hessian_operator_;
diff --git a/dune/gdt/operators/reconstruction/internal.hh b/dune/gdt/operators/reconstruction/internal.hh
index 850095c7c..787c2e16a 100644
--- a/dune/gdt/operators/reconstruction/internal.hh
+++ b/dune/gdt/operators/reconstruction/internal.hh
@@ -22,8 +22,6 @@
#include <dune/gdt/operators/interfaces.hh>
-#include "slopes.hh"
-
namespace Dune {
namespace GDT {
namespace internal {
diff --git a/dune/gdt/operators/reconstruction/linear_kinetic.hh b/dune/gdt/operators/reconstruction/linear_kinetic.hh
index 2fa385508..a2f55c91e 100644
--- a/dune/gdt/operators/reconstruction/linear_kinetic.hh
+++ b/dune/gdt/operators/reconstruction/linear_kinetic.hh
@@ -37,27 +37,22 @@ class LocalPointwiseLinearKineticReconstructionOperator : public XT::Grid::Eleme
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>;
+ static constexpr size_t stencil_size = 3;
+ using StencilType = XT::Common::FieldVector<size_t, stencil_size>;
+ using StencilsType = FieldVector<StencilType, dimDomain>;
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)
@@ -65,10 +60,7 @@ public:
, 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
@@ -84,22 +76,23 @@ public:
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);
+ analytical_flux_.calculate_reconstructed_fluxes(
+ stencils_[dd], boundary_dirs_[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];
+ stencils_[dd][1] = grid_view_.indexSet().index(entity);
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());
+ stencils_[dd][index] = size_t(-1);
+ boundary_dirs_[dd][index] = intersection.indexInInside() % 2;
} else if (intersection.neighbor()) { // inner and periodic intersections {
- stencils_[dd][index] = source_values_[grid_view_.indexSet().index(intersection.outside())];
+ stencils_[dd][index] = grid_view_.indexSet().index(intersection.outside());
} else if (!intersection.neighbor() && !intersection.boundary()) { // processor boundary
return false;
} else {
@@ -112,10 +105,10 @@ private:
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_;
+ XT::Common::FieldVector<XT::Common::FieldVector<size_t, stencil_size>, dimDomain> boundary_dirs_;
+
}; // class LocalPointwiseLinearKineticReconstructionOperator
// Does not reconstruct a full first-order DG function, but only stores the reconstructed values at the intersection
@@ -170,7 +163,7 @@ public:
// do reconstruction
auto local_reconstruction_operator =
LocalPointwiseLinearKineticReconstructionOperator<GV, AnalyticalFluxType, BoundaryValueType, LocalVectorType>(
- range, grid_view, analytical_flux_, source_values, boundary_values_, param);
+ range, grid_view, analytical_flux_, param);
auto walker = XT::Grid::Walker<GV>(grid_view);
walker.append(local_reconstruction_operator);
walker.walk(true);
diff --git a/dune/gdt/operators/reconstruction/slopes.hh b/dune/gdt/operators/reconstruction/slopes.hh
index a840a2aa7..4a3772eaa 100644
--- a/dune/gdt/operators/reconstruction/slopes.hh
+++ b/dune/gdt/operators/reconstruction/slopes.hh
@@ -31,8 +31,6 @@
#include <dune/gdt/momentmodels/basisfunctions/partial_moments.hh>
#include <dune/gdt/momentmodels/entropyflux.hh>
-#include "internal.hh"
-
namespace Dune {
namespace GDT {
@@ -213,22 +211,11 @@ 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.);
+ VectorType ret;
for (size_t ii = 0; ii < first_slope.size(); ++ii)
- if (std::signbit(first_slope[ii]) == std::signbit(second_slope[ii])) // check for equal sign
- ret[ii] = (std::abs(first_slope[ii]) < std::abs(second_slope[ii])) ? first_slope[ii] : second_slope[ii];
+ ret[ii] = XT::Common::minmod(first_slope[ii], second_slope[ii]);
return ret;
}
};
@@ -297,17 +284,14 @@ public:
VectorType slope_left_char, slope_right_char;
eigenvectors.apply_inverse_eigenvectors(dd, slope_left, slope_left_char);
eigenvectors.apply_inverse_eigenvectors(dd, slope_right, slope_right_char);
- const VectorType first_slope = MinmodType::minmod(slope_left_char, slope_right_char * 2.);
- const VectorType second_slope = MinmodType::minmod(slope_left_char * 2., slope_right_char);
- return maxmod(first_slope, second_slope);
+ return superbee(slope_left_char, slope_right_char);
}
- static VectorType maxmod(const VectorType& first_slope, const VectorType& second_slope)
+ static VectorType superbee(const VectorType& first_slope, const VectorType& second_slope)
{
VectorType ret(0.);
for (size_t ii = 0; ii < first_slope.size(); ++ii)
- if (std::signbit(first_slope[ii]) == std::signbit(second_slope[ii])) // check for equal sign
- ret[ii] = (std::abs(first_slope[ii]) < std::abs(second_slope[ii])) ? second_slope[ii] : first_slope[ii];
+ ret[ii] = superbee(first_slope[ii], second_slope[ii]);
return ret;
}
}; // class SuperbeeSlope
diff --git a/dune/gdt/test/entropic-coords-mn-discretization.hh b/dune/gdt/test/entropic-coords-mn-discretization.hh
index 220125d78..71e162969 100644
--- a/dune/gdt/test/entropic-coords-mn-discretization.hh
+++ b/dune/gdt/test/entropic-coords-mn-discretization.hh
@@ -28,6 +28,7 @@
#include <dune/gdt/momentmodels/hessianinverter.hh>
#include <dune/gdt/momentmodels/entropyflux_kineticcoords.hh>
#include <dune/gdt/momentmodels/entropyflux.hh>
+#include <dune/gdt/momentmodels/density_evaluations.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>
@@ -99,9 +100,11 @@ struct HyperbolicEntropicCoordsMnDiscretization
const auto& problem = *problem_ptr;
const auto initial_values_u = problem.initial_values();
const auto boundary_values_u = problem.boundary_values();
+ const auto boundary_distribution = problem.boundary_distribution();
using AnalyticalFluxType = typename ProblemType::FluxType;
- using EntropyFluxType = EntropyBasedFluxEntropyCoordsFunction<GV, MomentBasis, TestCaseType::reconstruction>;
+ constexpr SlopeType slope = TestCaseType::reconstruction ? SlopeType::minmod : SlopeType::no_slope;
+ using EntropyFluxType = EntropyBasedFluxEntropyCoordsFunction<GV, MomentBasis, slope>;
using OldEntropyFluxType = EntropyBasedFluxFunction<GV, MomentBasis>;
auto flux = problem.flux();
auto* entropy_flux = dynamic_cast<OldEntropyFluxType*>(flux.get());
@@ -123,18 +126,6 @@ struct HyperbolicEntropicCoordsMnDiscretization
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
@@ -153,12 +144,6 @@ struct HyperbolicEntropicCoordsMnDiscretization
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]);
}
@@ -166,7 +151,7 @@ struct HyperbolicEntropicCoordsMnDiscretization
// ******************** choose flux and rhs operator and timestepper ******************************************
using AdvectionOperatorType = AdvectionFvOperator<MatrixType, GV, dimRange>;
- using HessianInverterType = EntropicHessianInverter<MomentBasis, SpaceType, TestCaseType::reconstruction>;
+ using HessianInverterType = EntropicHessianInverter<MomentBasis, SpaceType, slope>;
#if 0
using ReconstructionOperatorType = PointwiseLinearReconstructionNoCharOperator<
GV,
@@ -180,14 +165,15 @@ struct HyperbolicEntropicCoordsMnDiscretization
#endif
using ReconstructionAdvectionOperatorType =
AdvectionWithPointwiseReconstructionOperator<AdvectionOperatorType, ReconstructionOperatorType>;
- using NonEntropicAdvectionOperatorType =
- std::conditional_t<TestCaseType::reconstruction, ReconstructionAdvectionOperatorType, AdvectionOperatorType>;
+ using NonEntropicAdvectionOperatorType = ReconstructionAdvectionOperatorType;
// using FvOperatorType = EntropicCoordinatesOperator<
// NonEntropicAdvectionOperatorType,
// HessianInverterType>;
using NonEntropicRhsOperatorType = GenericOperator<MatrixType, GV, dimRange>;
// using RhsOperatorType = EntropicCoordinatesOperator<NonEntropicRhsOperatorType, HessianInverterType>;
- using CombinedOperatorType = EntropicCoordinatesCombinedOperator<NonEntropicAdvectionOperatorType,
+ using DensityOperatorType = DensityEvaluator<MomentBasis, SpaceType, slope>;
+ using CombinedOperatorType = EntropicCoordinatesCombinedOperator<DensityOperatorType,
+ NonEntropicAdvectionOperatorType,
NonEntropicRhsOperatorType,
HessianInverterType>;
@@ -231,24 +217,38 @@ struct HyperbolicEntropicCoordsMnDiscretization
using LambdaType = typename BoundaryOperator::LambdaType;
using StateType = typename BoundaryOperator::StateType;
+ // calculate boundary kinetic fluxes
+ // apply density_operator first to store boundary_evaluations
+ const double min_acceptable_density = problem.psi_vac() / 10;
+ DensityOperatorType density_operator(*analytical_flux, fv_space, boundary_distribution, min_acceptable_density);
+ density_operator.apply(alpha.dofs().vector(), alpha.dofs().vector());
+
+ // store boundary 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 dd = intersection.indexInInside() / 2;
+ const auto boundary_flux = problem.kinetic_boundary_flux(x, intersection.centerUnitOuterNormal()[dd], dd);
+ 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]; });
+
+
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 TestCaseType::reconstruction
- ? boundary_kinetic_fluxes.evaluate(
- intersection.geometry().global(xx_in_reference_intersection_coordinates))
- : 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);
@@ -265,12 +265,8 @@ struct HyperbolicEntropicCoordsMnDiscretization
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);
-
+ HessianInverterType hessian_inverter(*analytical_flux, fv_space);
+ auto& non_entropic_advection_operator = reconstruction_advection_operator;
const auto u_iso = basis_functions->u_iso();
const auto basis_integrated = basis_functions->integrated();
@@ -278,14 +274,12 @@ struct HyperbolicEntropicCoordsMnDiscretization
const auto sigma_s = problem.sigma_s();
const auto Q = problem.Q();
auto rhs_func = [&](const auto& /*source*/,
- const auto& local_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 u = analytical_flux->get_u(fv_space.grid_view().indexSet().index(element));
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;
@@ -300,7 +294,7 @@ struct HyperbolicEntropicCoordsMnDiscretization
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);
+ density_operator, non_entropic_advection_operator, non_entropic_rhs_operator, hessian_inverter);
// ******************************** do the time steps ***********************************************************
// OperatorTimeStepperType timestepper_op(fv_operator, alpha, -1.0);
diff --git a/dune/gdt/test/hyperbolic__momentmodels__entropic_coords_mn.cc b/dune/gdt/test/hyperbolic__momentmodels__entropic_coords_mn.cc
index 0318d5c42..217629bbe 100644
--- a/dune/gdt/test/hyperbolic__momentmodels__entropic_coords_mn.cc
+++ b/dune/gdt/test/hyperbolic__momentmodels__entropic_coords_mn.cc
@@ -18,42 +18,32 @@ 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
+ // The kinetic scheme does not work for Legendre Mn in the SourceBeam test
+ // Dune::GDT::SourceBeamMnTestCase<Yasp1, Dune::GDT::LegendreMomentBasis<double, double, 7>, false, true>,
+ // Dune::GDT::SourceBeamMnTestCase<Yasp1, Dune::GDT::LegendreMomentBasis<double, double, 20>, true, true>
+ Dune::GDT::PlaneSourceMnTestCase<Yasp1, Dune::GDT::LegendreMomentBasis<double, double, 7>, false, true>,
+ Dune::GDT::PlaneSourceMnTestCase<Yasp1, Dune::GDT::LegendreMomentBasis<double, double, 7>, true, true>,
Dune::GDT::SourceBeamMnTestCase<Yasp1, Dune::GDT::HatFunctionMomentBasis<double, 1, double, 8, 1, 1>, false, true>,
- Dune::GDT::SourceBeamMnTestCase<Yasp1, Dune::GDT::HatFunctionMomentBasis<double, 1, double, 8, 1, 1>, true, true>
-// Dune::GDT::PlaneSourceMnTestCase<Yasp1, Dune::GDT::HatFunctionMomentBasis<double, 1, double, 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>
+ Dune::GDT::SourceBeamMnTestCase<Yasp1, Dune::GDT::HatFunctionMomentBasis<double, 1, double, 8, 1, 1>, true, true>,
+ Dune::GDT::PlaneSourceMnTestCase<Yasp1, Dune::GDT::HatFunctionMomentBasis<double, 1, double, 8, 1, 1>, false, true>,
+ Dune::GDT::PlaneSourceMnTestCase<Yasp1, Dune::GDT::HatFunctionMomentBasis<double, 1, double, 8, 1, 1>, true, true>,
+ Dune::GDT::SourceBeamMnTestCase<Yasp1, Dune::GDT::PartialMomentBasis<double, 1, double, 8, 1, 1>, false, true>,
+ Dune::GDT::SourceBeamMnTestCase<Yasp1, Dune::GDT::PartialMomentBasis<double, 1, double, 8, 1, 1>, true, true>,
+ Dune::GDT::PlaneSourceMnTestCase<Yasp1, Dune::GDT::PartialMomentBasis<double, 1, double, 8, 1, 1>, false, true>,
+ Dune::GDT::PlaneSourceMnTestCase<Yasp1, Dune::GDT::PartialMomentBasis<double, 1, double, 8, 1, 1>, true, 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
+ ,
+ Dune::GDT::
+ PointSourceMnTestCase<Yasp3, Dune::GDT::RealSphericalHarmonicsMomentBasis<double, double, 2, 3>, false, true>,
+ Dune::GDT::
+ PointSourceMnTestCase<Yasp3, Dune::GDT::RealSphericalHarmonicsMomentBasis<double, double, 2, 3>, true, true>,
+ // Dune::GDT::ShadowMnTestCase<Yasp3, Dune::GDT::HatFunctionMomentBasis<double, 3, double, 1, 1, 3>, true, true>,
+ Dune::GDT::CheckerboardMnTestCase<Yasp3, Dune::GDT::HatFunctionMomentBasis<double, 3, double, 0, 1, 3>, true, true>,
+ Dune::GDT::PointSourceMnTestCase<Yasp3, Dune::GDT::HatFunctionMomentBasis<double, 3, double, 0, 1, 3>, false, true>,
+ Dune::GDT::PointSourceMnTestCase<Yasp3, Dune::GDT::HatFunctionMomentBasis<double, 3, double, 0, 1, 3>, true, true>,
+ Dune::GDT::PointSourceMnTestCase<Yasp3, Dune::GDT::PartialMomentBasis<double, 3, double, 0, 1, 3>, false, true>,
+ Dune::GDT::PointSourceMnTestCase<Yasp3, Dune::GDT::PartialMomentBasis<double, 3, double, 0, 1, 3>, true, true>
+#endif // !DXT_DISABLE_LARGE_TESTS
>;
TYPED_TEST_CASE(HyperbolicEntropicCoordsMnTest, YaspGridTestCasesAll);
diff --git a/dune/gdt/test/momentmodels/kinetictransport/base.hh b/dune/gdt/test/momentmodels/kinetictransport/base.hh
index c092f43a9..e5f4165a8 100644
--- a/dune/gdt/test/momentmodels/kinetictransport/base.hh
+++ b/dune/gdt/test/momentmodels/kinetictransport/base.hh
@@ -57,6 +57,7 @@ public:
using DynamicFluxJacobianRangeType = typename FluxType::LocalFunctionType::DynamicJacobianRangeType;
using GenericScalarFunctionType = XT::Functions::GenericFunction<dimFlux, 1, 1, RangeFieldType>;
using ConstantScalarFunctionType = XT::Functions::ConstantFunction<dimFlux, 1, 1, RangeFieldType>;
+ using BoundaryDistributionType = std::function<std::function<RangeFieldType(const DomainType&)>(const DomainType&)>;
using BaseType::default_boundary_cfg;
using BaseType::default_grid_cfg;
@@ -173,6 +174,35 @@ public:
[=](const DomainType&, const XT::Common::Parameter&) { return value; });
} // ... boundary_values()
+ virtual BoundaryDistributionType boundary_distribution() const
+ {
+ return [this](const DomainType&) { return [this](const DomainType&) { return this->psi_vac_; }; };
+ }
+
+ RangeReturnType
+ kinetic_boundary_flux_from_quadrature(const DomainType& x, const RangeFieldType& n, const size_t dd) const
+ {
+ RangeReturnType ret(0.);
+ const auto boundary_density = boundary_distribution()(x);
+ const auto& quadratures = basis_functions_.quadratures();
+ for (size_t jj = 0; jj < quadratures.size(); ++jj) {
+ for (size_t ll = 0; ll < quadratures[jj].size(); ++ll) {
+ const auto v = quadratures[jj][ll].position();
+ const auto b = basis_functions_.evaluate(v, jj);
+ if (v[dd] * n < 0) {
+ const RangeFieldType psi = boundary_density(v);
+ ret += b * psi * v[dd] * quadratures[jj][ll].weight();
+ }
+ } // ll
+ } // jj
+ return ret;
+ }
+
+ virtual RangeReturnType kinetic_boundary_flux(const DomainType& x, const RangeFieldType& n, const size_t dd) const
+ {
+ return kinetic_boundary_flux_from_quadrature(x, n, dd);
+ }
+
virtual RangeFieldType CFL() const override
{
return 0.49;
diff --git a/dune/gdt/test/momentmodels/kinetictransport/sourcebeam.hh b/dune/gdt/test/momentmodels/kinetictransport/sourcebeam.hh
index adbeb01c9..93a22fd00 100644
--- a/dune/gdt/test/momentmodels/kinetictransport/sourcebeam.hh
+++ b/dune/gdt/test/momentmodels/kinetictransport/sourcebeam.hh
@@ -32,10 +32,12 @@ template <class E, class MomentBasisImp>
class SourceBeamPn : public KineticTransportEquationBase<E, MomentBasisImp>
{
using BaseType = KineticTransportEquationBase<E, MomentBasisImp>;
+ using ThisType = SourceBeamPn;
public:
using BaseType::dimDomain;
using BaseType::dimRange;
+ using typename BaseType::BoundaryDistributionType;
using typename BaseType::BoundaryValueType;
using typename BaseType::DomainFieldType;
using typename BaseType::DomainType;
@@ -92,6 +94,25 @@ public:
});
} // ... boundary_values()
+ virtual BoundaryDistributionType boundary_distribution() const override final
+ {
+ return [this](const DomainType& x) -> std::function<RangeFieldType(const DomainType&)> {
+ if (x[0] > 1.5)
+ return [this](const DomainType& /*v*/) { return this->psi_vac_; };
+ else
+ return [this](const DomainType& v) {
+ const RangeFieldType val = std::exp(-1e5 * std::pow(v[0] - 1., 2));
+ return (val > this->psi_vac_) ? val : this->psi_vac_;
+ };
+ };
+ }
+
+ virtual RangeReturnType
+ kinetic_boundary_flux(const DomainType& x, const RangeFieldType& n, const size_t dd) const override final
+ {
+ return helper<MomentBasis>::get_kinetic_boundary_flux(basis_functions_, psi_vac_, is_mn_model_, x, n, dd, *this);
+ }
+
RangeReturnType left_boundary_value() const
{
return helper<MomentBasis>::get_left_boundary_values(basis_functions_, psi_vac_, is_mn_model_);
@@ -133,11 +154,16 @@ protected:
struct helper_base
{
// returns the numerator g of the left boundary value (see create_boundary_values)
- static RangeFieldType numerator(const RangeFieldType v)
+ static RangeFieldType g(const RangeFieldType v)
{
return std::exp(-1e5 * (v - 1) * (v - 1));
}
+ static RangeFieldType dg(const RangeFieldType v)
+ {
+ return -2e5 * (v - 1) * g(v);
+ }
+
// returns the denominator <g> of the left boundary value (see create_boundary_values)
static RangeFieldType denominator()
{
@@ -155,7 +181,13 @@ protected:
// calculates integral from v_l to v_u of v*g
static RangeFieldType integral_2(RangeFieldType v_l, RangeFieldType v_u)
{
- return integral_1(v_l, v_u) - 1. / 2e5 * (numerator(v_u) - numerator(v_l));
+ return integral_1(v_l, v_u) - 1. / 2e5 * (g(v_u) - g(v_l));
+ }
+
+ // calculates integral from v_l to v_u of v^2*g
+ static RangeFieldType integral_3(RangeFieldType v_l, RangeFieldType v_u)
+ {
+ return 0.25e-10 * (dg(v_u) - dg(v_l)) + 2. * integral_2(v_l, v_u) + (0.5e-5 - 1.) * integral_1(v_l, v_u);
}
};
@@ -164,7 +196,7 @@ protected:
struct helper : public helper_base
{
using helper_base::denominator;
- using helper_base::numerator;
+ using helper_base::g;
static DynamicRangeType get_left_boundary_values(const MomentBasisImp& basis_functions,
const RangeFieldType& psi_vac,
@@ -182,7 +214,7 @@ protected:
for (const auto& quad_point : quadrature) {
const auto& v = quad_point.position()[0];
auto summand = basis_functions.evaluate(v, ii);
- summand *= numerator(v) * quad_point.weight();
+ summand *= g(v) * quad_point.weight();
ret += summand;
}
}
@@ -191,6 +223,17 @@ protected:
ret += basis_functions.integrated() * psi_vac;
return ret;
}
+
+ static DynamicRangeType get_kinetic_boundary_flux(const MomentBasisImp& /*basis_functions*/,
+ const RangeFieldType& /*psi_vac*/,
+ const bool /*is_mn_model*/,
+ const DomainType& x,
+ const RangeFieldType& n,
+ const size_t dd,
+ const ThisType& problem)
+ {
+ return problem.kinetic_boundary_flux_from_quadrature(x, n, dd);
+ }
};
template <class anything>
@@ -198,8 +241,10 @@ protected:
: public helper_base
{
using helper_base::denominator;
+ using helper_base::g;
using helper_base::integral_1;
- using helper_base::numerator;
+ using helper_base::integral_2;
+ using helper_base::integral_3;
static DynamicRangeType get_left_boundary_values(const MomentBasisImp& basis_functions,
const RangeFieldType psi_vac,
@@ -211,19 +256,55 @@ protected:
const auto vn = triangulation[nn];
if (nn < dimRange - 1) {
const auto vnp = triangulation[nn + 1];
- ret[nn] += 1. / ((vn - vnp) * denominator())
- * ((1 - vnp) * integral_1(vn, vnp) - 1. / 2e5 * (numerator(vnp) - numerator(vn)));
+ ret[nn] += 1. / ((vn - vnp) * denominator()) * (integral_2(vn, vnp) - vnp * integral_1(vn, vnp));
}
if (nn > 0) {
const auto vnm = triangulation[nn - 1];
- ret[nn] += 1. / ((vn - vnm) * denominator())
- * ((1 - vnm) * integral_1(vnm, vn) - 1. / 2e5 * (numerator(vn) - numerator(vnm)));
+ ret[nn] += 1. / ((vn - vnm) * denominator()) * (integral_2(vnm, vn) - vnm * integral_1(vnm, vn));
}
}
// add small vacuum concentration to move away from realizable boundary
ret += basis_functions.integrated() * psi_vac;
return ret;
- }
+ } // ... get_left_boundary_values(...)
+
+ static DynamicRangeType get_kinetic_boundary_flux(const MomentBasisImp& basis_functions,
+ const RangeFieldType& /*psi_vac*/,
+ const bool is_mn_model,
+ const DomainType& x,
+ const RangeFieldType& n,
+ const size_t dd,
+ const ThisType& problem)
+ {
+ if (!is_mn_model)
+ DUNE_THROW(Dune::NotImplemented, "Only implemented for mn");
+ if (x < 1.5) {
+ DynamicRangeType ret(dimRange, 0.);
+ const auto& triangulation = basis_functions.triangulation();
+ for (size_t nn = 0; nn < dimRange; ++nn) {
+ const auto vn = triangulation[nn];
+ if (nn < dimRange - 1) {
+ const auto vnp = triangulation[nn + 1];
+ if (vnp > 0.) {
+ const auto left_limit = vn > 0. ? vn : 0.;
+ ret[nn] +=
+ 1. / ((vn - vnp) * denominator()) * (integral_3(left_limit, vnp) - vnp * integral_2(left_limit, vnp));
+ } // if (vnp > 0.)
+ } // if (nn < dimRange -1)
+ if (vn > 0.) {
+ if (nn > 0) {
+ const auto vnm = triangulation[nn - 1];
+ const auto left_limit = vnm > 0. ? vnm : 0.;
+ ret[nn] +=
+ 1. / ((vn - vnm) * denominator()) * (integral_2(left_limit, vn) - vnm * integral_1(left_limit, vn));
+ } // if (nn > 0)
+ } // if (vn > 0.)
+ } // nn
+ return ret;
+ } else {
+ return problem.kinetic_boundary_flux_from_quadrature(x, n, dd);
+ }
+ } // ... get_kinetic_boundary_flux(...)
};
template <class anything>
@@ -232,6 +313,7 @@ protected:
using helper_base::denominator;
using helper_base::integral_1;
using helper_base::integral_2;
+ using helper_base::integral_3;
static DynamicRangeType get_left_boundary_values(const MomentBasisImp& basis_functions,
const RangeFieldType psi_vac,
@@ -247,6 +329,32 @@ protected:
ret += basis_functions.integrated() * psi_vac;
return ret;
}
+
+ static DynamicRangeType get_kinetic_boundary_flux(const MomentBasisImp& basis_functions,
+ const RangeFieldType& /*psi_vac*/,
+ const bool is_mn_model,
+ const DomainType& x,
+ const RangeFieldType& n,
+ const size_t dd,
+ const ThisType& problem)
+ {
+ if (!is_mn_model)
+ DUNE_THROW(Dune::NotImplemented, "Only implemented for mn");
+ if (x < 1.5) {
+ const auto& triangulation = basis_functions.triangulation();
+ DynamicRangeType ret(dimRange, 0.);
+ for (size_t ii = 0; ii < dimRange / 2; ++ii) {
+ if (triangulation[ii + 1] > 0.) {
+ const auto left_limit = triangulation[ii] > 0. ? triangulation[ii] : 0.;
+ ret[2 * ii] = integral_2(left_limit, triangulation[ii + 1]) / denominator();
+ ret[2 * ii + 1] = integral_3(left_limit, triangulation[ii + 1]) / denominator();
+ }
+ } // ii
+ return ret;
+ } else {
+ return problem.kinetic_boundary_flux_from_quadrature(x, n, dd);
+ }
+ } // ... get_kinetic_boundary_flux(...)
};
using BaseType::basis_functions_;
diff --git a/dune/gdt/test/momentmodels/kinetictransport/testcases.hh b/dune/gdt/test/momentmodels/kinetictransport/testcases.hh
index bc5b7bd3f..11b99cfb2 100644
--- a/dune/gdt/test/momentmodels/kinetictransport/testcases.hh
+++ b/dune/gdt/test/momentmodels/kinetictransport/testcases.hh
@@ -238,7 +238,7 @@ struct SourceBeamPnTestCase
// SourceBeam Mn
-template <class MomentBasisImp, bool reconstruct, bool kinetic_scheme>
+template <class MomentBasisImp, bool reconstruct, bool kinetic_scheme = false>
struct SourceBeamMnExpectedResults;
template <bool reconstruct, bool kinetic_scheme>
@@ -262,14 +262,14 @@ struct SourceBeamMnExpectedResults<HatFunctionMomentBasis<double, 1, double, 8,
template <bool reconstruct>
struct SourceBeamMnExpectedResults<HatFunctionMomentBasis<double, 1, double, 8, 1, 1>, reconstruct, true>
{
- static constexpr double l1norm = reconstruct ? 370.93823722076462 : 367.56820428572496;
- static constexpr double l2norm = reconstruct ? 235.97872826934559 : 235.22990290508343;
- static constexpr double linfnorm = reconstruct ? 210.82343000666447 : 209.02263326452186;
+ static constexpr double l1norm = reconstruct ? 371.54588397717055 : 367.97988291905477;
+ static constexpr double l2norm = reconstruct ? 236.4476851910448 : 235.54814675091959;
+ static constexpr double linfnorm = reconstruct ? 210.63369526083264 : 208.81107020771216;
static constexpr double tol = 1e-9;
};
-template <bool reconstruct, bool kinetic_scheme>
-struct SourceBeamMnExpectedResults<PartialMomentBasis<double, 1, double, 8, 1, 1>, reconstruct, kinetic_scheme>
+template <bool reconstruct>
+struct SourceBeamMnExpectedResults<PartialMomentBasis<double, 1, double, 8, 1, 1>, reconstruct, false>
{
static constexpr double l1norm = reconstruct ? 0.33140398337368543 : 0.3314039833756291;
static constexpr double l2norm = reconstruct ? 0.45583354074069732 : 0.44484887610818585;
@@ -277,6 +277,15 @@ struct SourceBeamMnExpectedResults<PartialMomentBasis<double, 1, double, 8, 1, 1
static constexpr double tol = 1e-9;
};
+template <bool reconstruct>
+struct SourceBeamMnExpectedResults<PartialMomentBasis<double, 1, double, 8, 1, 1>, reconstruct, true>
+{
+ static constexpr double l1norm = reconstruct ? 254.20216502516391 : 270.74268779687191;
+ static constexpr double l2norm = reconstruct ? 187.86036790841933 : 202.76054800096165;
+ static constexpr double linfnorm = reconstruct ? 265.10790627160509 : 260.82089045524185;
+ static constexpr double tol = 1e-9;
+};
+
template <class GridImp, class MomentBasisImp, bool reconstruct, bool kinetic_scheme = false>
struct SourceBeamMnTestCase : public SourceBeamPnTestCase<GridImp, MomentBasisImp, reconstruct>
{
@@ -334,11 +343,17 @@ struct PlaneSourcePnTestCase : SourceBeamPnTestCase<GridImp, MomentBasisImp, rec
// PlaneSource Mn
-template <class MomentBasisImp, bool reconstruct>
-struct PlaneSourceMnExpectedResults;
+template <class MomentBasisImp, bool reconstruct, bool kinetic_scheme = false>
+struct PlaneSourceMnExpectedResults
+{
+ static constexpr double l1norm = 0.;
+ static constexpr double l2norm = 0.;
+ static constexpr double linfnorm = 0.;
+ static constexpr double tol = 0.;
+};
template <bool reconstruct>
-struct PlaneSourceMnExpectedResults<LegendreMomentBasis<double, double, 7>, reconstruct>
+struct PlaneSourceMnExpectedResults<LegendreMomentBasis<double, double, 7>, reconstruct, false>
{
static constexpr double l1norm = reconstruct ? 2.0000000240000007 : 2.0000000240000029;
static constexpr double l2norm = reconstruct ? 2.785411193059216 : 2.746101358507282;
@@ -347,7 +362,16 @@ struct PlaneSourceMnExpectedResults<LegendreMomentBasis<double, double, 7>, reco
};
template <bool reconstruct>
-struct PlaneSourceMnExpectedResults<HatFunctionMomentBasis<double, 1, double, 8, 1, 1>, reconstruct>
+struct PlaneSourceMnExpectedResults<LegendreMomentBasis<double, double, 7>, reconstruct, true>
+{
+ static constexpr double l1norm = reconstruct ? 33.830651291425575 : 31.119878976551046;
+ static constexpr double l2norm = reconstruct ? 24.726893737746675 : 23.385570207485049;
+ static constexpr double linfnorm = reconstruct ? 19.113827924512311 : 19.113827924512311;
+ static constexpr double tol = 1e-7;
+};
+
+template <bool reconstruct>
+struct PlaneSourceMnExpectedResults<HatFunctionMomentBasis<double, 1, double, 8, 1, 1>, reconstruct, false>
{
static constexpr double l1norm = 2.0000000239315696;
static constexpr double l2norm = reconstruct ? 2.7966600752714887 : 2.7457411547488615;
@@ -356,7 +380,16 @@ struct PlaneSourceMnExpectedResults<HatFunctionMomentBasis<double, 1, double, 8,
};
template <bool reconstruct>
-struct PlaneSourceMnExpectedResults<PartialMomentBasis<double, 1, double, 8, 1, 1>, reconstruct>
+struct PlaneSourceMnExpectedResults<HatFunctionMomentBasis<double, 1, double, 8, 1, 1>, reconstruct, true>
+{
+ static constexpr double l1norm = reconstruct ? 268.42786311776274 : 246.7429359648828;
+ static constexpr double l2norm = reconstruct ? 197.1506642519208 : 186.09403264481648;
+ static constexpr double linfnorm = reconstruct ? 152.91062339609854 : 152.91062339609854;
+ static constexpr double tol = 1e-9;
+};
+
+template <bool reconstruct>
+struct PlaneSourceMnExpectedResults<PartialMomentBasis<double, 1, double, 8, 1, 1>, reconstruct, false>
{
static constexpr double l1norm = reconstruct ? 2.0000000239999913 : 2.0000000239999904;
static constexpr double l2norm = reconstruct ? 2.8215879031834015 : 2.7633864171098814;
@@ -364,7 +397,17 @@ struct PlaneSourceMnExpectedResults<PartialMomentBasis<double, 1, double, 8, 1,
static constexpr double tol = 1e-9;
};
-template <class GridImp, class MomentBasisImp, bool reconstruct>
+template <bool reconstruct>
+struct PlaneSourceMnExpectedResults<PartialMomentBasis<double, 1, double, 8, 1, 1>, reconstruct, true>
+{
+ static constexpr double l1norm = reconstruct ? 144.19157186249112 : 135.86834797834712;
+ static constexpr double l2norm = reconstruct ? 104.28938402311856 : 100.2359224660796;
+ static constexpr double linfnorm = reconstruct ? 100.43185554232102 : 97.933985765677008;
+ static constexpr double tol = 1e-9;
+};
+
+
+template <class GridImp, class MomentBasisImp, bool reconstruct, bool kinetic_scheme = false>
struct PlaneSourceMnTestCase : SourceBeamMnTestCase<GridImp, MomentBasisImp, reconstruct>
{
using BaseType = SourceBeamMnTestCase<GridImp, MomentBasisImp, reconstruct>;
@@ -374,7 +417,7 @@ struct PlaneSourceMnTestCase : SourceBeamMnTestCase<GridImp, MomentBasisImp, rec
using ProblemType = PlaneSourceMn<GridViewType, MomentBasisImp>;
static constexpr RangeFieldType t_end = 0.25;
static constexpr bool reconstruction = reconstruct;
- using ExpectedResultsType = PlaneSourceMnExpectedResults<MomentBasisImp, reconstruction>;
+ using ExpectedResultsType = PlaneSourceMnExpectedResults<MomentBasisImp, reconstruction, kinetic_scheme>;
using RealizabilityLimiterChooserType =
RealizabilityLimiterChooser<GridViewType, MomentBasisImp, typename ProblemType::FluxType, DiscreteFunctionType>;
};
@@ -522,11 +565,17 @@ struct ShadowPnTestCase : SourceBeamPnTestCase<GridImp, MomentBasisImp, reconstr
// PointSourceMn
-template <class MomentBasisImp, bool reconstruct>
-struct PointSourceMnExpectedResults;
+template <class MomentBasisImp, bool reconstruct, bool kinetic_scheme = false>
+struct PointSourceMnExpectedResults
+{
+ static constexpr double l1norm = 0.;
+ static constexpr double l2norm = 0.;
+ static constexpr double linfnorm = 0.;
+ static constexpr double tol = 0.;
+};
template <bool reconstruct>
-struct PointSourceMnExpectedResults<RealSphericalHarmonicsMomentBasis<double, double, 2, 3>, reconstruct>
+struct PointSourceMnExpectedResults<RealSphericalHarmonicsMomentBasis<double, double, 2, 3>, reconstruct, false>
{
static constexpr double l1norm = reconstruct ? 1.0000013830443908 : 1.0000013830442143;
static constexpr double l2norm = reconstruct ? 2.6901467570598112 : 2.684314243798307;
@@ -535,7 +584,16 @@ struct PointSourceMnExpectedResults<RealSphericalHarmonicsMomentBasis<double, do
};
template <bool reconstruct>
-struct PointSourceMnExpectedResults<HatFunctionMomentBasis<double, 3, double, 0, 1, 3>, reconstruct>
+struct PointSourceMnExpectedResults<RealSphericalHarmonicsMomentBasis<double, double, 2, 3>, reconstruct, true>
+{
+ static constexpr double l1norm = reconstruct ? 1674.9008041695579 : 1585.7044225325101;
+ static constexpr double l2norm = reconstruct ? 619.41343145125302 : 589.93299566257235;
+ static constexpr double linfnorm = reconstruct ? 264.16080528868997 : 266.5453598444696;
+ static constexpr double tol = 1e-9;
+};
+
+template <bool reconstruct>
+struct PointSourceMnExpectedResults<HatFunctionMomentBasis<double, 3, double, 0, 1, 3>, reconstruct, false>
{
static constexpr double l1norm = reconstruct ? 1.0000000829624791 : 1.0000000829622864;
static constexpr double l2norm = reconstruct ? 2.694751941188763 : 2.6892684619955305;
@@ -551,7 +609,16 @@ struct PointSourceMnExpectedResults<HatFunctionMomentBasis<double, 3, double, 0,
};
template <bool reconstruct>
-struct PointSourceMnExpectedResults<PartialMomentBasis<double, 3, double, 0, 1, 3, 1>, reconstruct>
+struct PointSourceMnExpectedResults<HatFunctionMomentBasis<double, 3, double, 0, 1, 3>, reconstruct, true>
+{
+ static constexpr double l1norm = reconstruct ? 743.43592672927934 : 683.47651349520368;
+ static constexpr double l2norm = reconstruct ? 279.11700895978396 : 258.45009663177717;
+ static constexpr double linfnorm = reconstruct ? 125.7102296804655 : 125.71014355518233;
+ static constexpr double tol = 1e-9;
+};
+
+template <bool reconstruct>
+struct PointSourceMnExpectedResults<PartialMomentBasis<double, 3, double, 0, 1, 3, 1>, reconstruct, false>
{
static constexpr double l1norm = reconstruct ? 1.0000000829624787 : 1.0000000829623072;
static constexpr double l2norm = reconstruct ? 2.6983516853120966 : 2.6881937835020211;
@@ -559,16 +626,25 @@ struct PointSourceMnExpectedResults<PartialMomentBasis<double, 3, double, 0, 1,
static constexpr double tol = 1e-9;
};
-template <class GridImp, class MomentBasisImp, bool reconstruct>
-struct PointSourceMnTestCase : SourceBeamMnTestCase<GridImp, MomentBasisImp, reconstruct>
+template <bool reconstruct>
+struct PointSourceMnExpectedResults<PartialMomentBasis<double, 3, double, 0, 1, 3, 1>, reconstruct, true>
{
- using BaseType = SourceBeamMnTestCase<GridImp, MomentBasisImp, reconstruct>;
+ static constexpr double l1norm = reconstruct ? 1167.5985275432627 : 1126.5174600848904;
+ static constexpr double l2norm = reconstruct ? 428.15220455351266 : 415.19451310103005;
+ static constexpr double linfnorm = reconstruct ? 188.26590907577059 : 191.48910050886076;
+ static constexpr double tol = 1e-9;
+};
+
+template <class GridImp, class MomentBasisImp, bool reconstruct, bool kinetic_scheme = false>
+struct PointSourceMnTestCase : SourceBeamMnTestCase<GridImp, MomentBasisImp, reconstruct, kinetic_scheme>
+{
+ using BaseType = SourceBeamMnTestCase<GridImp, MomentBasisImp, reconstruct, kinetic_scheme>;
using typename BaseType::GridViewType;
using ProblemType = PointSourceMn<GridViewType, MomentBasisImp>;
using typename BaseType::RangeFieldType;
static constexpr RangeFieldType t_end = 0.1;
static constexpr bool reconstruction = reconstruct;
- using ExpectedResultsType = PointSourceMnExpectedResults<MomentBasisImp, reconstruction>;
+ using ExpectedResultsType = PointSourceMnExpectedResults<MomentBasisImp, reconstruction, kinetic_scheme>;
using RealizabilityLimiterChooserType = RealizabilityLimiterChooser<GridViewType,
MomentBasisImp,
typename ProblemType::FluxType,
@@ -577,11 +653,11 @@ struct PointSourceMnTestCase : SourceBeamMnTestCase<GridImp, MomentBasisImp, rec
// CheckerboardMn
-template <class MomentBasisImp, bool reconstruct>
+template <class MomentBasisImp, bool reconstruct, bool kinetic_scheme = false>
struct CheckerboardMnExpectedResults;
template <bool reconstruct>
-struct CheckerboardMnExpectedResults<RealSphericalHarmonicsMomentBasis<double, double, 2, 3>, reconstruct>
+struct CheckerboardMnExpectedResults<RealSphericalHarmonicsMomentBasis<double, double, 2, 3>, reconstruct, false>
{
static constexpr double l1norm = reconstruct ? 0. : 0.35404509573284748;
static constexpr double l2norm = reconstruct ? 0. : 0.32922954029850499;
@@ -589,16 +665,43 @@ struct CheckerboardMnExpectedResults<RealSphericalHarmonicsMomentBasis<double, d
static constexpr double tol = 1e-9;
};
-template <class GridImp, class MomentBasisImp, bool reconstruct>
-struct CheckerboardMnTestCase : SourceBeamMnTestCase<GridImp, MomentBasisImp, reconstruct>
+template <bool reconstruct>
+struct CheckerboardMnExpectedResults<RealSphericalHarmonicsMomentBasis<double, double, 2, 3>, reconstruct, true>
{
- using BaseType = SourceBeamMnTestCase<GridImp, MomentBasisImp, reconstruct>;
+ static constexpr double l1norm = reconstruct ? 0. : 0.;
+ static constexpr double l2norm = reconstruct ? 0. : 0.;
+ static constexpr double linfnorm = reconstruct ? 0. : 0.;
+ static constexpr double tol = 1e-9;
+};
+
+template <bool reconstruct>
+struct CheckerboardMnExpectedResults<HatFunctionMomentBasis<double, 3, double, 0, 1, 3>, reconstruct, false>
+{
+ static constexpr double l1norm = reconstruct ? 0. : 0.;
+ static constexpr double l2norm = reconstruct ? 0. : 0.;
+ static constexpr double linfnorm = reconstruct ? 0. : 0.;
+ static constexpr double tol = 1e-9;
+};
+
+template <bool reconstruct>
+struct CheckerboardMnExpectedResults<HatFunctionMomentBasis<double, 3, double, 0, 1, 3>, reconstruct, true>
+{
+ static constexpr double l1norm = reconstruct ? 44760.517221844952 : 0.;
+ static constexpr double l2norm = reconstruct ? 2491.7035345029099 : 0.;
+ static constexpr double linfnorm = reconstruct ? 542.02861997383934 : 0.;
+ static constexpr double tol = 1e-9;
+};
+
+template <class GridImp, class MomentBasisImp, bool reconstruct, bool kinetic_scheme = false>
+struct CheckerboardMnTestCase : SourceBeamMnTestCase<GridImp, MomentBasisImp, reconstruct, kinetic_scheme>
+{
+ using BaseType = SourceBeamMnTestCase<GridImp, MomentBasisImp, reconstruct, kinetic_scheme>;
using typename BaseType::GridViewType;
using ProblemType = CheckerboardMn<GridViewType, MomentBasisImp>;
using typename BaseType::RangeFieldType;
static constexpr RangeFieldType t_end = 0.1;
static constexpr bool reconstruction = reconstruct;
- using ExpectedResultsType = CheckerboardMnExpectedResults<MomentBasisImp, reconstruction>;
+ using ExpectedResultsType = CheckerboardMnExpectedResults<MomentBasisImp, reconstruction, kinetic_scheme>;
using RealizabilityLimiterChooserType = RealizabilityLimiterChooser<GridViewType,
MomentBasisImp,
typename ProblemType::FluxType,
--
GitLab