diff --git a/dune/gdt/local/numerical-fluxes/kinetic.hh b/dune/gdt/local/numerical-fluxes/kinetic.hh
index 925c30b608b7f73b09acab4257b30df970701b79..5783fecb35ba1de40ffc0451a3a22fcd0399d0a3 100644
--- a/dune/gdt/local/numerical-fluxes/kinetic.hh
+++ b/dune/gdt/local/numerical-fluxes/kinetic.hh
@@ -48,17 +48,17 @@ class NumericalKineticFlux
public:
using typename BaseType::FluxType;
- using typename BaseType::FunctionType;
using typename BaseType::LocalIntersectionCoords;
using typename BaseType::PhysicalDomainType;
using typename BaseType::StateType;
+ using typename BaseType::XIndependentFluxType;
NumericalKineticFlux(const FluxType& flx, const MomentBasis& basis)
: BaseType(flx)
, basis_(basis)
{}
- NumericalKineticFlux(const FunctionType& flx, const MomentBasis& basis)
+ NumericalKineticFlux(const XIndependentFluxType& flx, const MomentBasis& basis)
: BaseType(flx)
, basis_(basis)
{}
@@ -72,20 +72,20 @@ public:
using BaseType::apply;
using BaseType::flux;
+ using BaseType::intersection;
- StateType apply(const I& intersection,
- const LocalIntersectionCoords& /*x*/,
+ StateType apply(const LocalIntersectionCoords& /*x*/,
const StateType& u,
const StateType& v,
const PhysicalDomainType& n,
const XT::Common::Parameter& /*param*/ = {}) const override final
{
// find direction of unit outer normal (we assume an axis-aligned cube grid)
- size_t direction = intersection.indexInInside() / 2;
+ size_t direction = intersection().indexInInside() / 2;
if (dynamic_cast<const EntropyFluxType*>(&flux()) != nullptr) {
return dynamic_cast<const EntropyFluxType*>(&flux())->evaluate_kinetic_flux(
- intersection.inside(),
- intersection.neighbor() ? intersection.outside() : intersection.inside(),
+ intersection().inside(),
+ intersection().neighbor() ? intersection().outside() : intersection().inside(),
u,
v,
n,
diff --git a/dune/gdt/momentmodels/basisfunctions/hatfunctions.hh b/dune/gdt/momentmodels/basisfunctions/hatfunctions.hh
index af900d860191ec59e7c32ea0979d865456097b75..680f7749e15b0351191bc0c689f71e1382916ab1 100644
--- a/dune/gdt/momentmodels/basisfunctions/hatfunctions.hh
+++ b/dune/gdt/momentmodels/basisfunctions/hatfunctions.hh
@@ -63,6 +63,13 @@ public:
return "hatfunctions";
}
+ static size_t default_quad_order()
+ {
+ return 15;
+ }
+
+ using BaseType::default_quad_refinements;
+
HatFunctionMomentBasis(const QuadraturesType& quadratures)
: BaseType(quadratures)
, triangulation_(BaseType::create_1d_triangulation(dimRange - 1))
@@ -70,7 +77,8 @@ public:
BaseType::initialize_base_values();
}
- HatFunctionMomentBasis(const size_t quad_order = 15, const size_t DXTC_DEBUG_ONLY(quad_refinements) = 0)
+ HatFunctionMomentBasis(const size_t quad_order = default_quad_order(),
+ const size_t DXTC_DEBUG_ONLY(quad_refinements) = default_quad_refinements())
: BaseType(BaseType::gauss_lobatto_quadratures(dimRange - 1, quad_order))
, triangulation_(BaseType::create_1d_triangulation(dimRange - 1))
{
@@ -355,6 +363,13 @@ public:
using BaseType::barycentre_rule;
+ static size_t default_quad_order()
+ {
+ return refinements == 0 ? 15 : 9;
+ }
+
+ using BaseType::default_quad_refinements;
+
HatFunctionMomentBasis(const QuadraturesType& quadratures)
: BaseType(refinements, quadratures)
{
@@ -371,7 +386,8 @@ public:
BaseType::initialize_base_values();
}
- HatFunctionMomentBasis(const size_t quad_order = (refinements == 0 ? 15 : 9), const size_t quad_refinements = 0)
+ HatFunctionMomentBasis(const size_t quad_order = default_quad_order(),
+ const size_t quad_refinements = default_quad_refinements())
: BaseType(refinements)
{
const QuadratureRule<RangeFieldType, 2> reference_quadrature_rule =
diff --git a/dune/gdt/momentmodels/basisfunctions/interface.hh b/dune/gdt/momentmodels/basisfunctions/interface.hh
index af132d836433ddf4dcc3f2160013bfed989c6ddf..6b76300201b0ffc952de51cf4edb9c3367eb5116 100644
--- a/dune/gdt/momentmodels/basisfunctions/interface.hh
+++ b/dune/gdt/momentmodels/basisfunctions/interface.hh
@@ -168,6 +168,17 @@ public:
using MergedQuadratureIterator =
typename XT::Data::MergedQuadrature<RangeFieldType, dimDomain>::MergedQuadratureIterator;
+ static size_t default_quad_order()
+ {
+ DUNE_THROW(Dune::NotImplemented, "Please overwrite this function in derived classes!");
+ return 0;
+ }
+
+ static size_t default_quad_refinements()
+ {
+ return 0;
+ }
+
MomentBasisInterface(const QuadraturesType& quadratures = QuadraturesType())
: quadratures_(quadratures)
{}
diff --git a/dune/gdt/momentmodels/basisfunctions/legendre.hh b/dune/gdt/momentmodels/basisfunctions/legendre.hh
index 93162bf940b6db4b58c55e75ddb57aa17185d987..8aad82d2b0c85d92ac2bc0b4353c566dafc1878c 100644
--- a/dune/gdt/momentmodels/basisfunctions/legendre.hh
+++ b/dune/gdt/momentmodels/basisfunctions/legendre.hh
@@ -42,7 +42,18 @@ public:
BaseType::initialize_base_values();
}
- LegendreMomentBasis(const size_t quad_order = 197, const size_t quad_refinements = 1)
+ static size_t default_quad_order()
+ {
+ return 2 * order + 40;
+ }
+
+ static size_t default_quad_refinements()
+ {
+ return 1;
+ }
+
+ LegendreMomentBasis(const size_t quad_order = default_quad_order(),
+ const size_t quad_refinements = default_quad_refinements())
: BaseType(BaseType::gauss_lobatto_quadratures(std::pow(2, quad_refinements), quad_order))
{
BaseType::initialize_base_values();
diff --git a/dune/gdt/momentmodels/basisfunctions/partial_moments.hh b/dune/gdt/momentmodels/basisfunctions/partial_moments.hh
index b38dc638d340b3d6c24e844ba8339796ab72b15f..a8531bf02c97ed0896a2ac3d716e7ccde55515a4 100644
--- a/dune/gdt/momentmodels/basisfunctions/partial_moments.hh
+++ b/dune/gdt/momentmodels/basisfunctions/partial_moments.hh
@@ -70,6 +70,13 @@ public:
return "pcw";
}
+ static size_t default_quad_order()
+ {
+ return 15;
+ }
+
+ using BaseType::default_quad_refinements;
+
PartialMomentBasis(const QuadraturesType& quadratures)
: BaseType(quadratures)
, triangulation_(BaseType::create_1d_triangulation(num_intervals))
@@ -77,7 +84,8 @@ public:
BaseType::initialize_base_values();
}
- PartialMomentBasis(const size_t quad_order = 15, const size_t DXTC_DEBUG_ONLY(quad_refinements) = 0)
+ PartialMomentBasis(const size_t quad_order = default_quad_order(),
+ const size_t DXTC_DEBUG_ONLY(quad_refinements) = default_quad_refinements())
: BaseType(BaseType::gauss_lobatto_quadratures(num_intervals, quad_order))
, triangulation_(BaseType::create_1d_triangulation(num_intervals))
{
@@ -406,6 +414,13 @@ public:
using BaseType::barycentre_rule;
+ static size_t default_quad_order()
+ {
+ return refinements == 0 ? 15 : 9;
+ }
+
+ using BaseType::default_quad_refinements;
+
PartialMomentBasis(const QuadraturesType& quadratures)
: BaseType(refinements, quadratures)
{
@@ -421,7 +436,8 @@ public:
BaseType::initialize_base_values();
}
- PartialMomentBasis(const size_t quad_order = (refinements == 0 ? 15 : 9), const size_t quad_refinements = 0)
+ PartialMomentBasis(const size_t quad_order = default_quad_order(),
+ const size_t quad_refinements = default_quad_refinements())
: BaseType(refinements)
{
const QuadratureRule<RangeFieldType, 2> reference_quadrature_rule =
diff --git a/dune/gdt/momentmodels/basisfunctions/spherical_harmonics.hh b/dune/gdt/momentmodels/basisfunctions/spherical_harmonics.hh
index c13c036dfe9d244ef46ff9f781384499eccf219a..ad590e5d58161c03f710dc2c16aecbe4d3129319 100644
--- a/dune/gdt/momentmodels/basisfunctions/spherical_harmonics.hh
+++ b/dune/gdt/momentmodels/basisfunctions/spherical_harmonics.hh
@@ -50,14 +50,21 @@ public:
using typename BaseType::StringifierType;
static_assert(order <= std::numeric_limits<int>::max(), "");
+ static size_t default_quad_order()
+ {
+ return 2 * order + 8;
+ }
+
+ using BaseType::default_quad_refinements;
+
SphericalHarmonicsMomentBasis(const QuadraturesType& quadratures)
: BaseType(quadratures)
{
BaseType::initialize_base_values();
}
- SphericalHarmonicsMomentBasis(const size_t quad_order = 2 * order + 8,
- const size_t DXTC_DEBUG_ONLY(quad_refinements) = 0)
+ SphericalHarmonicsMomentBasis(const size_t quad_order = default_quad_order(),
+ const size_t DXTC_DEBUG_ONLY(quad_refinements) = default_quad_refinements())
: BaseType(XT::Data::OctantQuadratures<DomainFieldType>::get(quad_order))
{
assert(quad_refinements == 0 && "Refinement of the quadrature intervals not implemented for this basis!");
@@ -289,14 +296,21 @@ public:
using typename BaseType::RangeType;
using typename BaseType::StringifierType;
+ static size_t default_quad_order()
+ {
+ return 2 * order + 8;
+ }
+
+ using BaseType::default_quad_refinements;
+
RealSphericalHarmonicsMomentBasis(const QuadraturesType& quadratures)
: BaseType(quadratures)
{
BaseType::initialize_base_values();
}
- RealSphericalHarmonicsMomentBasis(const size_t quad_order = 2 * order + 8,
- const size_t DXTC_DEBUG_ONLY(quad_refinements) = 0)
+ RealSphericalHarmonicsMomentBasis(const size_t quad_order = default_quad_order(),
+ const size_t DXTC_DEBUG_ONLY(quad_refinements) = default_quad_refinements())
: BaseType(XT::Data::OctantQuadratures<DomainFieldType>::get(quad_order))
{
assert(quad_refinements == 0 && "Refinement of the quadrature intervals not implemented for this basis!");
diff --git a/dune/gdt/momentmodels/entropyflux.hh b/dune/gdt/momentmodels/entropyflux.hh
index 52b43473c60be76b3e49761ae7188d2cdbdc42ff..2710961bb8a9dd85135a79a4d267a39f5eee3ace 100644
--- a/dune/gdt/momentmodels/entropyflux.hh
+++ b/dune/gdt/momentmodels/entropyflux.hh
@@ -157,11 +157,33 @@ public:
const size_t k_max = 1000,
const RangeFieldType epsilon = std::pow(2, -52))
: index_set_(grid_view.indexSet())
+ , use_thread_cache_(true)
+ , use_entity_cache_(true)
, entity_caches_(index_set_.size(0), LocalCacheType(cache_size))
, mutexes_(index_set_.size(0))
, implementation_(basis_functions, tau, epsilon_gamma, chi, xi, r_sequence, k_0, k_max, epsilon)
{}
+ void enable_thread_cache()
+ {
+ use_thread_cache_ = true;
+ }
+
+ void disable_thread_cache()
+ {
+ use_thread_cache_ = false;
+ }
+
+ void enable_entity_cache()
+ {
+ use_entity_cache_ = true;
+ }
+
+ void disable_entity_cache()
+ {
+ use_entity_cache_ = false;
+ }
+
static const constexpr bool available = true;
class Localfunction : public LocalFunctionType
@@ -175,11 +197,15 @@ public:
Localfunction(const IndexSetType& index_set,
std::vector<LocalCacheType>& entity_caches,
+ const bool use_thread_cache,
+ const bool use_entity_cache,
std::vector<std::mutex>& mutexes,
const ImplementationType& implementation)
: index_set_(index_set)
, thread_cache_(cache_size)
, entity_caches_(entity_caches)
+ , use_thread_cache_(use_thread_cache)
+ , use_entity_cache_(use_entity_cache)
, mutexes_(mutexes)
, implementation_(implementation)
{}
@@ -209,7 +235,7 @@ public:
// calculate (inf-norm) distance to isotropic moment with same density
RangeFieldType distance = (u - u_iso_scaled).infinity_norm();
VectorType alpha_start = XT::Common::convert_to<VectorType>(alpha_iso + alpha_iso_prime * std::log(density));
- if (!XT::Common::FloatCmp::eq(distance, 0.)) {
+ if (!XT::Common::FloatCmp::eq(distance, 0.) && use_entity_cache_) {
// calculate distance to closest moment in entity_cache
const auto entity_cache_dist_and_it = entity_cache_->find_closest(u);
const auto& entity_cache_dist = entity_cache_dist_and_it.first;
@@ -217,7 +243,7 @@ public:
distance = entity_cache_dist;
alpha_start = entity_cache_dist_and_it.second->second;
}
- if (!XT::Common::FloatCmp::eq(distance, 0.)) {
+ if (!XT::Common::FloatCmp::eq(distance, 0.) && use_thread_cache_) {
// calculate distance to closest moment in thread_cache
const auto thread_cache_dist_and_it = thread_cache_.find_closest(u);
const auto& thread_cache_dist = thread_cache_dist_and_it.first;
@@ -232,8 +258,10 @@ public:
return std::make_unique<AlphaReturnType>(std::make_pair(alpha_start, std::make_pair(u, 0.)));
} else {
auto ret = implementation_.get_alpha(u, alpha_start, regularize);
- entity_cache_->insert(ret->second.first, ret->first);
- thread_cache_.insert(ret->second.first, ret->first);
+ if (use_entity_cache_)
+ entity_cache_->insert(ret->second.first, ret->first);
+ if (use_thread_cache_)
+ thread_cache_.insert(ret->second.first, ret->first);
return std::move(ret);
}
}
@@ -259,6 +287,8 @@ public:
const IndexSetType& index_set_;
mutable LocalCacheType thread_cache_;
std::vector<LocalCacheType>& entity_caches_;
+ const bool use_thread_cache_;
+ const bool use_entity_cache_;
std::vector<std::mutex>& mutexes_;
const ImplementationType& implementation_;
LocalCacheType* entity_cache_;
@@ -272,12 +302,14 @@ public:
virtual std::unique_ptr<LocalFunctionType> local_function() const override final
{
- return std::make_unique<Localfunction>(index_set_, entity_caches_, mutexes_, implementation_);
+ return std::make_unique<Localfunction>(
+ index_set_, entity_caches_, use_thread_cache_, use_entity_cache_, mutexes_, implementation_);
}
virtual std::unique_ptr<Localfunction> derived_local_function() const
{
- return std::make_unique<Localfunction>(index_set_, entity_caches_, mutexes_, implementation_);
+ return std::make_unique<Localfunction>(
+ index_set_, entity_caches_, use_thread_cache_, use_entity_cache_, mutexes_, implementation_);
}
StateType evaluate_kinetic_flux(const E& inside_entity,
@@ -303,6 +335,8 @@ public:
private:
const IndexSetType& index_set_;
+ bool use_thread_cache_;
+ bool use_entity_cache_;
mutable std::vector<LocalCacheType> entity_caches_;
mutable std::vector<std::mutex> mutexes_;
ImplementationType implementation_;
diff --git a/dune/gdt/momentmodels/entropyflux_implementations.hh b/dune/gdt/momentmodels/entropyflux_implementations.hh
index 5f25b1fcd39f8e7af47a0b11e3a20f2e8d9ec306..ed141053fb423640626b0a71da87045edafe23c2 100644
--- a/dune/gdt/momentmodels/entropyflux_implementations.hh
+++ b/dune/gdt/momentmodels/entropyflux_implementations.hh
@@ -463,7 +463,6 @@ protected:
1);
#else
const size_t num_quad_points = quad_points_.size();
- std::fill(scalar_products.begin(), scalar_products.end(), 0.);
for (size_t ll = 0; ll < num_quad_points; ++ll) {
const auto* basis_ll = &(M.get_entry_ref(ll, 0.));
scalar_products[ll] = std::inner_product(beta_in.begin(), beta_in.end(), basis_ll, 0.);
@@ -743,7 +742,7 @@ public:
// regularize u
v = u_prime;
if (r > 0) {
- beta_in = get_isotropic_alpha(u);
+ beta_in = get_isotropic_alpha(v);
VectorType r_times_u_iso = u_iso;
r_times_u_iso *= r;
v *= 1 - r;
@@ -942,7 +941,7 @@ public:
// regularize u
v = u_prime;
if (r > 0) {
- alpha_k = get_isotropic_alpha(u);
+ alpha_k = get_isotropic_alpha(v);
VectorType r_times_u_iso = u_iso;
r_times_u_iso *= r;
v *= 1 - r;
@@ -1276,7 +1275,7 @@ public:
// regularize u
*v = *u_prime;
if (r > 0.) {
- *beta_in = *get_isotropic_alpha(u);
+ *beta_in = *get_isotropic_alpha(*v);
// calculate v = (1-r) u + r u_iso
// use beta_out as storage for u_iso_in * r
*v *= (1 - r);
@@ -1848,6 +1847,14 @@ public:
return ret;
}
+ VectorType get_isotropic_alpha(const VectorType& u) const
+ {
+ DomainType u_domain;
+ for (size_t ii = 0; ii < basis_dimDomain; ++ii)
+ u_domain[ii] = u.get_entry(ii);
+ return get_isotropic_alpha(u_domain);
+ }
+
virtual RangeReturnType evaluate(const DomainType& u,
const XT::Common::Parameter& /*param*/ = {}) const override final
{
@@ -1985,7 +1992,7 @@ public:
// regularize u
v = u_prime;
if (r > 0) {
- alpha_k = get_isotropic_alpha(u);
+ alpha_k = get_isotropic_alpha(v);
tmp_vec = u_iso;
tmp_vec *= r;
v *= 1 - r;
@@ -2615,7 +2622,7 @@ public:
// regularize u
v = u_prime;
if (r > 0) {
- alpha_k = get_isotropic_alpha(u);
+ alpha_k = get_isotropic_alpha(v);
VectorType r_times_u_iso(u_iso);
r_times_u_iso *= r;
v *= 1 - r;
@@ -3023,16 +3030,15 @@ public:
using QuadraturePointsType = FieldVector<std::vector<RangeFieldType>, num_intervals>;
using QuadratureWeightsType = FieldVector<std::vector<RangeFieldType>, num_intervals>;
- explicit EntropyBasedFluxImplementation(
- const MomentBasis& basis_functions,
- const RangeFieldType tau = 1e-9,
- const RangeFieldType epsilon_gamma = 0.01,
- const RangeFieldType chi = 0.5,
- const RangeFieldType xi = 1e-3,
- const std::vector<RangeFieldType> r_sequence = {0, 1e-8, 1e-6, 1e-4, 1e-3, 1e-2, 5e-2, 0.1, 0.5, 1},
- const size_t k_0 = 500,
- const size_t k_max = 1000,
- const RangeFieldType epsilon = std::pow(2, -52))
+ explicit EntropyBasedFluxImplementation(const MomentBasis& basis_functions,
+ const RangeFieldType tau,
+ const RangeFieldType epsilon_gamma,
+ const RangeFieldType chi,
+ const RangeFieldType xi,
+ const std::vector<RangeFieldType> r_sequence,
+ const size_t k_0,
+ const size_t k_max,
+ const RangeFieldType epsilon)
: basis_functions_(basis_functions)
, grid_points_(basis_functions_.triangulation())
, tau_(tau)
@@ -3186,7 +3192,7 @@ public:
// regularize u
v = u_prime;
if (r > 0) {
- alpha_k = get_isotropic_alpha(u);
+ alpha_k = get_isotropic_alpha(v);
VectorType r_times_u_iso(u_iso);
r_times_u_iso *= r;
v *= 1 - r;
diff --git a/dune/gdt/momentmodels/entropysolver.hh b/dune/gdt/momentmodels/entropysolver.hh
index bb1cbf4cb7b581d86261188a2c9525da77e3695a..c43f48c3204047ace559c0e3d7180b6df95dfd94 100644
--- a/dune/gdt/momentmodels/entropysolver.hh
+++ b/dune/gdt/momentmodels/entropysolver.hh
@@ -88,8 +88,6 @@ public:
for (size_t ii = 0; ii < dimRange; ++ii)
u[ii] = local_source_->dofs().get_entry(ii);
const auto& basis_functions = analytical_flux_.basis_functions();
- thread_local auto vector_indices = source_.space().mapper().global_indices(entity);
- source_.space().mapper().global_indices(entity, vector_indices);
basis_functions.ensure_min_density(u, min_acceptable_density_);
local_flux_->bind(entity);
const auto alpha_ret = local_flux_->get_alpha(u, true);
diff --git a/dune/gdt/operators/reconstruction/internal.hh b/dune/gdt/operators/reconstruction/internal.hh
index d65f93d949448ce5d9603cc55b8232a89c34a826..f2453d752604ff9d5ad22567df5a82d042f64801 100644
--- a/dune/gdt/operators/reconstruction/internal.hh
+++ b/dune/gdt/operators/reconstruction/internal.hh
@@ -194,9 +194,9 @@ public:
const XT::Common::Parameter& param) override final
{
local_flux_->bind(entity);
- local_flux_->jacobian(x_local, u, *jacobian_, param);
- for (size_t dd = 0; dd < dimDomain; ++dd) {
- try {
+ try {
+ local_flux_->jacobian(x_local, u, *jacobian_, param);
+ for (size_t dd = 0; dd < dimDomain; ++dd) {
if (false) {
;
#if HAVE_MKL || HAVE_LAPACKE
@@ -249,14 +249,17 @@ public:
if (info || eigenvectors_rcond_ < 1e-5)
DUNE_THROW(Dune::MathError, "Eigenvector condition too high!");
#endif
- } catch (const Dune::MathError&) {
+ } // dd
+ } catch (const Dune::MathError&) {
+ for (size_t dd = 0; dd < dimDomain; ++dd) {
// use scalar limiters, i.e. eigenvectors matrix is eye-matrix.
XT::LA::eye_matrix((*eigenvectors_)[dd]);
std::fill(eigenvalues_[dd].begin(), eigenvalues_[dd].end(), 1.);
(*QR_)[dd] = (*eigenvectors_)[dd];
XT::LA::qr((*QR_)[dd], tau_[dd], permutations_[dd]);
- }
- } // dd
+ } // dd
+ }
+
// we do not need jacobian_ anymore if the flux is affine, so save memory
if (flux_is_affine_)
jacobian_ = nullptr;
@@ -343,11 +346,11 @@ public:
{
local_flux_->bind(entity);
const FluxDomainType nonblocked_u = u.operator FluxDomainType();
- local_flux_->jacobian(x_local, nonblocked_u, *nonblocked_jacobian_, param);
- *jacobian_ = *nonblocked_jacobian_;
- for (size_t dd = 0; dd < dimDomain; ++dd) {
- for (size_t jj = 0; jj < num_blocks; ++jj) {
- try {
+ try {
+ local_flux_->jacobian(x_local, nonblocked_u, *nonblocked_jacobian_, param);
+ *jacobian_ = *nonblocked_jacobian_;
+ for (size_t dd = 0; dd < dimDomain; ++dd) {
+ for (size_t jj = 0; jj < num_blocks; ++jj) {
if (block_size == 2) {
const auto& jac = (*jacobian_)[dd].block(jj);
const auto trace = jac[0][0] + jac[1][1];
@@ -414,15 +417,20 @@ public:
if (info || eigenvectors_rcond_ < 1e-5)
DUNE_THROW(Dune::MathError, "Eigenvector condition too high!");
#endif
- } catch (const Dune::MathError&) {
+ } // jj
+ } // dd
+ } catch (const Dune::MathError&) {
+ for (size_t dd = 0; dd < dimDomain; ++dd) {
+ for (size_t jj = 0; jj < num_blocks; ++jj) {
// use scalar limiters, i.e. eigenvectors matrix is eye-matrix.
XT::LA::eye_matrix(eigenvectors_[dd].block(jj));
std::fill(eigenvalues_[dd][jj].begin(), eigenvalues_[dd][jj].end(), 1.);
QR_[dd].block(jj) = eigenvectors_[dd].block(jj);
XT::LA::qr(QR_[dd].block(jj), tau_[dd].block(jj), permutations_[dd].block(jj));
- }
- } // jj
- } // dd
+ } // jj
+ } // dd
+ }
+
if (flux_is_affine_) {
jacobian_ = nullptr;
nonblocked_jacobian_ = nullptr;
diff --git a/dune/gdt/operators/reconstruction/linear.hh b/dune/gdt/operators/reconstruction/linear.hh
index b350a297cd2bc4f3af544a3e72181ce919fce139..375214d16ad36f5fcb2236124d5350547ae3aea6 100644
--- a/dune/gdt/operators/reconstruction/linear.hh
+++ b/dune/gdt/operators/reconstruction/linear.hh
@@ -14,7 +14,6 @@
#include <dune/geometry/quadraturerules.hh>
#include <dune/xt/common/fvector.hh>
-#include <dune/xt/common/parallel/threadstorage.hh>
#include <dune/xt/common/parameter.hh>
#include <dune/xt/grid/walker.hh>
@@ -396,6 +395,9 @@ private:
}; // class LinearReconstructionOperator<...>
+// Does not reconstruct a full first-order DG function, but only stores the reconstructed values at the intersection
+// centers. This avoids the interpolation in this operator and the evaluation of the reconstructed function in the
+// finite volume operator which are both quite expensive for large dimRange.
template <class AnalyticalFluxImp,
class BoundaryValueImp,
class GV,
diff --git a/dune/gdt/operators/reconstruction/slopes.hh b/dune/gdt/operators/reconstruction/slopes.hh
index b1533e323be68f3781fba6ac6aaab30b471aabe7..3c94056799254ae3f03e60a78431f6cd7da5b5b6 100644
--- a/dune/gdt/operators/reconstruction/slopes.hh
+++ b/dune/gdt/operators/reconstruction/slopes.hh
@@ -1011,40 +1011,32 @@ public:
private:
void setup_linear_program() const
{
- // It should be possible to reuse the same ClpSimplex class over and over again, only changing the relevant columns.
- // However, in rare cases the class seems to get corrupted, always claiming the LP is infeasible, even if it is not.
- // Creating a new LP from time to time seems to fix this problem.
- thread_local int counter;
- ++counter;
- if (!lp_ || !(counter % 100)) {
- counter = 0;
- // We start with creating a model with dimRange rows and num_quad_points+1 columns */
- assert(basis_values_->size() + dimRange < std::numeric_limits<int>::max());
- size_t num_cols = basis_values_->size() + dimRange; /* variables are x_1, ..., x_{num_quad_points}, theta_1,
- ..., theta_{dimRange} */
- lp_ = std::make_unique<ClpSimplex>(false);
- auto& lp = *lp_;
- // set number of rows
- lp.resize(num_rows, 0);
-
- // set columns for quadrature points
- assert(basis_values_->size() == num_cols - dimRange);
- static auto row_indices = create_row_indices();
- for (size_t ii = 0; ii < num_cols - dimRange; ++ii) {
- // First argument: number of elements in column
- // Second/Third argument: indices/values of column entries
- // Fourth/Fifth argument: lower/upper column bound, i.e. lower/upper bound for x_i. As all x_i should be
- // positive, set to 0/inf, which is the default.
- // Sixth argument: Prefactor in objective for x_i, this is 0 for all x_i, which is also the default;
- lp.addColumn(static_cast<int>(num_rows), row_indices.data(), &((*basis_values_)[ii][0]));
- }
+ // We start with creating a model with dimRange rows and num_quad_points+1 columns */
+ assert(basis_values_->size() + dimRange < std::numeric_limits<int>::max());
+ size_t num_cols = basis_values_->size() + dimRange; /* variables are x_1, ..., x_{num_quad_points}, theta_1,
+ ..., theta_{dimRange} */
+ lp_ = std::make_unique<ClpSimplex>(false);
+ auto& lp = *lp_;
+ // set number of rows
+ lp.resize(num_rows, 0);
+
+ // set columns for quadrature points
+ assert(basis_values_->size() == num_cols - dimRange);
+ static auto row_indices = create_row_indices();
+ for (size_t ii = 0; ii < num_cols - dimRange; ++ii) {
+ // First argument: number of elements in column
+ // Second/Third argument: indices/values of column entries
+ // Fourth/Fifth argument: lower/upper column bound, i.e. lower/upper bound for x_i. As all x_i should be
+ // positive, set to 0/inf, which is the default.
+ // Sixth argument: Prefactor in objective for x_i, this is 0 for all x_i, which is also the default;
+ lp.addColumn(static_cast<int>(num_rows), row_indices.data(), &((*basis_values_)[ii][0]));
+ }
- // add theta columns (set to random values, will be set correctly in solve_linear_program)
- // The bounds for theta should be [0,1]. Also sets the prefactor in the objective to 1 for the thetas.
- for (size_t ii = 0; ii < dimRange; ++ii)
- lp.addColumn(static_cast<int>(num_rows), row_indices.data(), &((*basis_values_)[0][0]), 0., 1., 1.);
- lp.setLogLevel(0);
- } // if (!lp_)
+ // add theta columns (set to random values, will be set correctly in solve_linear_program)
+ // The bounds for theta should be [0,1]. Also sets the prefactor in the objective to 1 for the thetas.
+ for (size_t ii = 0; ii < dimRange; ++ii)
+ lp.addColumn(static_cast<int>(num_rows), row_indices.data(), &((*basis_values_)[0][0]), 0., 1., 1.);
+ lp.setLogLevel(0);
} // void setup_linear_program()
VectorType solve_linear_program(const VectorType& u_char,
diff --git a/dune/gdt/test/mn-discretization.hh b/dune/gdt/test/mn-discretization.hh
index f0cae0222e4a5f884e8a6895474a79744cc15056..55c670d42a254cf8080defac9edb9aeae082088d 100644
--- a/dune/gdt/test/mn-discretization.hh
+++ b/dune/gdt/test/mn-discretization.hh
@@ -39,14 +39,15 @@ template <class TestCaseType>
struct HyperbolicMnDiscretization
{
// returns: (l1norm, l2norm, linfnorm, MPI rank)
- static std::pair<Dune::FieldVector<double, 3>, int>
- run(size_t num_save_steps = 1,
- size_t num_output_steps = 0,
- size_t quad_order = TestCaseType::RealizabilityLimiterChooserType::quad_order,
- std::string grid_size = "",
- std::string overlap_size = "2",
- double t_end = TestCaseType::t_end,
- std::string filename = "")
+ static std::pair<Dune::FieldVector<double, 3>, int> run(size_t num_save_steps = 1,
+ size_t num_output_steps = 0,
+ size_t quad_order = size_t(-1),
+ size_t quad_refinements = size_t(-1),
+ size_t grid_size = size_t(-1),
+ size_t overlap_size = 2,
+ double t_end = 0.,
+ std::string filename = "",
+ bool disable_thread_cache = false)
{
using namespace Dune;
using namespace Dune::GDT;
@@ -69,9 +70,9 @@ struct HyperbolicMnDiscretization
//******************* create grid and FV space ***************************************
auto grid_config = ProblemType::default_grid_cfg();
- if (!grid_size.empty())
- grid_config["num_elements"] = grid_size;
- grid_config["overlap_size"] = overlap_size;
+ if (grid_size != size_t(-1))
+ grid_config["num_elements"] = XT::Common::to_string(grid_size);
+ grid_config["overlap_size"] = XT::Common::to_string(overlap_size);
const auto grid_ptr =
Dune::XT::Grid::CubeGridProviderFactory<GridType>::create(grid_config, MPIHelper::getCommunicator()).grid_ptr();
assert(grid_ptr->comm().size() == 1 || grid_ptr->overlapSize(0) > 0);
@@ -80,7 +81,9 @@ struct HyperbolicMnDiscretization
const AdvectionSourceSpaceType advection_source_space(grid_view, 1);
//******************* create EquationType object ***************************************
- std::shared_ptr<const MomentBasis> basis_functions = std::make_shared<const MomentBasis>(quad_order);
+ std::shared_ptr<const MomentBasis> basis_functions = std::make_shared<const MomentBasis>(
+ quad_order == size_t(-1) ? MomentBasis::default_quad_order() : quad_order,
+ quad_refinements == size_t(-1) ? MomentBasis::default_quad_refinements() : quad_refinements);
const std::unique_ptr<ProblemType> problem_ptr =
XT::Common::make_unique<ProblemType>(*basis_functions, grid_view, grid_config);
const auto& problem = *problem_ptr;
@@ -88,7 +91,12 @@ struct HyperbolicMnDiscretization
const auto boundary_values = problem.boundary_values();
using AnalyticalFluxType = typename ProblemType::FluxType;
using EntropyFluxType = typename ProblemType::ActualFluxType;
- const auto analytical_flux = problem.flux();
+ auto analytical_flux = problem.flux();
+ // for Legendre polynomials and real spherical harmonics, the results are sensitive to the initial guess in the
+ // Newton algorithm. If the thread cache is enabled, the guess is different dependent on how many threads we are
+ // using, so for the tests we disable this cache to get reproducible results.
+ if (disable_thread_cache)
+ dynamic_cast<EntropyFluxType*>(analytical_flux.get())->disable_thread_cache();
const RangeFieldType CFL = problem.CFL();
// ***************** project initial values to discrete function *********************
@@ -153,18 +161,23 @@ struct HyperbolicMnDiscretization
constexpr double epsilon = 1e-11;
auto slope = TestCaseType::RealizabilityLimiterChooserType::template make_slope<EigenvectorWrapperType>(
- *dynamic_cast<const EntropyFluxType*>(analytical_flux.get()), *basis_functions, epsilon);
+ *dynamic_cast<EntropyFluxType*>(analytical_flux.get()), *basis_functions, epsilon);
ReconstructionOperatorType reconstruction_operator(*analytical_flux, *boundary_values, fv_space, *slope, false);
ReconstructionAdvectionOperatorType reconstruction_advection_operator(advection_operator, reconstruction_operator);
- filename += "_" + ProblemType::static_id();
- filename += "_grid_" + grid_size;
+ if (XT::Common::is_zero(t_end))
+ t_end = problem.t_end();
+
+ if (!filename.empty())
+ filename += "_";
+ filename += ProblemType::static_id();
+ filename += "_grid_" + grid_config["num_elements"];
filename += "_tend_" + XT::Common::to_string(t_end);
filename += "_quad_" + XT::Common::to_string(quad_order);
filename += TestCaseType::reconstruction ? "_ord2" : "_ord1";
filename += "_" + basis_functions->mn_name();
- EntropySolverType entropy_solver(*(dynamic_cast<const EntropyFluxType*>(analytical_flux.get())),
+ EntropySolverType entropy_solver(*(dynamic_cast<EntropyFluxType*>(analytical_flux.get())),
fv_space,
problem.psi_vac() * basis_functions->unit_ball_volume() / 10,
filename);
@@ -226,7 +239,17 @@ struct HyperbolicMnTest
{
void run()
{
- auto norms = HyperbolicMnDiscretization<TestCaseType>::run().first;
+ auto norms = HyperbolicMnDiscretization<TestCaseType>::run(
+ 1,
+ 0,
+ TestCaseType::RealizabilityLimiterChooserType::quad_order,
+ TestCaseType::RealizabilityLimiterChooserType::quad_refinements,
+ size_t(-1),
+ 2,
+ TestCaseType::t_end,
+ "test",
+ Dune::GDT::is_full_moment_basis<typename TestCaseType::MomentBasis>::value)
+ .first;
const double l1norm = norms[0];
const double l2norm = norms[1];
const double linfnorm = norms[2];
diff --git a/dune/gdt/test/momentmodels/kinetictransport/testcases.hh b/dune/gdt/test/momentmodels/kinetictransport/testcases.hh
index 9b84bd06a65ca22c6b6346d66e5e9ed8dddc4232..12ba0aa3e2f003c67b63bbba3be4f48d26cb856d 100644
--- a/dune/gdt/test/momentmodels/kinetictransport/testcases.hh
+++ b/dune/gdt/test/momentmodels/kinetictransport/testcases.hh
@@ -43,6 +43,7 @@ struct RealizabilityLimiterChooser<GV,
using MomentBasis = LegendreMomentBasis<double, double, order>;
using EntropyFluxType = EntropyBasedFluxFunction<GV, MomentBasis>;
static constexpr size_t quad_order = 54;
+ static constexpr size_t quad_refinements = 1;
template <class EigenVectorWrapperType>
static std::unique_ptr<LpConvexhullRealizabilityLimitedSlope<GV, MomentBasis, EigenVectorWrapperType>>
@@ -66,6 +67,7 @@ struct RealizabilityLimiterChooser<GV,
using MomentBasis = HatFunctionMomentBasis<double, 1, double, dimRange, 1, 1>;
using EntropyFluxType = EntropyBasedFluxFunction<GV, MomentBasis>;
static constexpr size_t quad_order = 15;
+ static constexpr size_t quad_refinements = 0;
#if HAVE_CLP && USE_LP_POSITIVITY_LIMITER
template <class EigenVectorWrapperType>
@@ -95,6 +97,7 @@ struct RealizabilityLimiterChooser<GV,
using MomentBasis = PartialMomentBasis<double, 1, double, dimRange, 1, 1>;
using EntropyFluxType = EntropyBasedFluxFunction<GV, MomentBasis>;
static constexpr size_t quad_order = 15;
+ static constexpr size_t quad_refinements = 0;
template <class EigenVectorWrapperType>
static std::unique_ptr<Dg1dRealizabilityLimitedSlope<GV, double, dimRange, EigenVectorWrapperType>>
@@ -114,7 +117,8 @@ struct RealizabilityLimiterChooser<GV,
{
using MomentBasis = RealSphericalHarmonicsMomentBasis<double, double, order, 3>;
using EntropyFluxType = EntropyBasedFluxFunction<GV, MomentBasis>;
- static constexpr size_t quad_order = 2 * order + 6;
+ static constexpr size_t quad_order = 2 * order + 8;
+ static constexpr size_t quad_refinements = 0;
template <class EigenVectorWrapperType>
static std::unique_ptr<LpConvexhullRealizabilityLimitedSlope<GV, MomentBasis, EigenVectorWrapperType>>
@@ -135,7 +139,8 @@ struct RealizabilityLimiterChooser<GV,
using MomentBasis = HatFunctionMomentBasis<double, 3, double, refinements, 1, 3>;
using EntropyFluxType = EntropyBasedFluxFunction<GV, MomentBasis>;
static constexpr size_t dimRange = MomentBasis::dimRange;
- static constexpr size_t quad_order = 18; // fekete rule number 7
+ static constexpr size_t quad_order = refinements == 0 ? 18 /*fekete rule number 7*/ : 9 /*fekete rule number 3*/;
+ static constexpr size_t quad_refinements = 0;
#if HAVE_CLP && USE_LP_POSITIVITY_LIMITER
template <class EigenVectorWrapperType>
@@ -165,7 +170,8 @@ struct RealizabilityLimiterChooser<GV,
{
using MomentBasis = PartialMomentBasis<double, 3, double, refinements, 1, 3>;
using EntropyFluxType = EntropyBasedFluxFunction<GV, MomentBasis>;
- static constexpr size_t quad_order = 9; // fekete rule number 3
+ static constexpr size_t quad_order = refinements == 0 ? 18 /*fekete rule number 7*/ : 9 /*fekete rule number 3*/;
+ static constexpr size_t quad_refinements = 0;
template <class EigenVectorWrapperType>
static std::unique_ptr<DgConvexHullRealizabilityLimitedSlope<GV, MomentBasis, EigenVectorWrapperType>>
@@ -239,9 +245,9 @@ template <bool reconstruct>
struct SourceBeamMnExpectedResults<LegendreMomentBasis<double, double, 7>, reconstruct>
{
static constexpr double l1norm = reconstruct ? 0.28535354296013105 : 0.28535354295945792;
- static constexpr double l2norm = reconstruct ? 0.37115083996060916 : 0.36265752973701221;
+ static constexpr double l2norm = reconstruct ? 0.37115145999473981 : 0.36265752973701221;
static constexpr double linfnorm = reconstruct ? 0.78506610334488358 : 0.78315544039143314;
- static constexpr double tol = 1e-9;
+ static constexpr double tol = 1e-5;
};
template <bool reconstruct>
@@ -326,9 +332,9 @@ template <bool reconstruct>
struct PlaneSourceMnExpectedResults<LegendreMomentBasis<double, double, 7>, reconstruct>
{
static constexpr double l1norm = reconstruct ? 2.0000000240000007 : 2.0000000240000029;
- static constexpr double l2norm = reconstruct ? 2.7878232892168211 : 2.746101358507282;
- static constexpr double linfnorm = reconstruct ? 4.9086815941726396 : 5.327698357914608;
- static constexpr double tol = 1e-9;
+ static constexpr double l2norm = reconstruct ? 2.785411193059216 : 2.746101358507282;
+ static constexpr double linfnorm = reconstruct ? 4.9069101475812698 : 5.327698357914608;
+ static constexpr double tol = 1e-7;
};
template <bool reconstruct>
@@ -513,9 +519,9 @@ struct PointSourceMnExpectedResults;
template <bool reconstruct>
struct PointSourceMnExpectedResults<RealSphericalHarmonicsMomentBasis<double, double, 2, 3>, reconstruct>
{
- static constexpr double l1norm = reconstruct ? 1.0000016440585395 : 1.0000016440583621;
- static constexpr double l2norm = reconstruct ? 2.6896356393271734 : 2.6837985426830504;
- static constexpr double linfnorm = reconstruct ? 10.368884563879304 : 10.375079241728713;
+ static constexpr double l1norm = reconstruct ? 1.0000013830443908 : 1.0000013830442143;
+ static constexpr double l2norm = reconstruct ? 2.6901467570598112 : 2.684314243798307;
+ static constexpr double linfnorm = reconstruct ? 10.371048798431969 : 10.377307670780343;
static constexpr double tol = 1e-9;
};
@@ -538,9 +544,9 @@ struct PointSourceMnExpectedResults<HatFunctionMomentBasis<double, 3, double, 0,
template <bool reconstruct>
struct PointSourceMnExpectedResults<PartialMomentBasis<double, 3, double, 0, 1, 3, 1>, reconstruct>
{
- static constexpr double l1norm = reconstruct ? 1.0000127558028837 : 1.0000127558027134;
- static constexpr double l2norm = reconstruct ? 2.6983746461758269 : 2.6882147627660715;
- static constexpr double linfnorm = reconstruct ? 10.391230462003152 : 10.394186618150204;
+ static constexpr double l1norm = reconstruct ? 1.0000000829624787 : 1.0000000829623072;
+ static constexpr double l2norm = reconstruct ? 2.6983516853120966 : 2.6881937835020211;
+ static constexpr double linfnorm = reconstruct ? 10.391142640527102 : 10.394108065213185;
static constexpr double tol = 1e-9;
};
@@ -568,9 +574,9 @@ struct CheckerboardMnExpectedResults;
template <bool reconstruct>
struct CheckerboardMnExpectedResults<RealSphericalHarmonicsMomentBasis<double, double, 2, 3>, reconstruct>
{
- static constexpr double l1norm = reconstruct ? 0. : 0.35404367426148931;
- static constexpr double l2norm = reconstruct ? 0. : 0.32921585327065972;
- static constexpr double linfnorm = reconstruct ? 0. : 0.32895499053937993;
+ static constexpr double l1norm = reconstruct ? 0. : 0.35404440392013337;
+ static constexpr double l2norm = reconstruct ? 0. : 0.32922954029850499;
+ static constexpr double linfnorm = reconstruct ? 0. : 0.32896894056609421;
static constexpr double tol = 1e-9;
};
@@ -598,9 +604,9 @@ struct ShadowMnExpectedResults;
template <bool reconstruct>
struct ShadowMnExpectedResults<RealSphericalHarmonicsMomentBasis<double, double, 2, 3>, reconstruct>
{
- static constexpr double l1norm = reconstruct ? 0. : 0.59247415692031025;
- static constexpr double l2norm = reconstruct ? 0. : 0.097642936357459895;
- static constexpr double linfnorm = reconstruct ? 0. : 0.016480937983211638;
+ static constexpr double l1norm = reconstruct ? 0. : 0.59248402251960053;
+ static constexpr double l2norm = reconstruct ? 0. : 0.097644561106262767;
+ static constexpr double linfnorm = reconstruct ? 0. : 0.016480889201743513;
static constexpr double tol = 1e-9;
};
diff --git a/dune/gdt/test/pn-discretization.hh b/dune/gdt/test/pn-discretization.hh
index 6f3ca032250c171af1098075258a47e88b7b2cab..29f144ec478dde96044d748d2a8886d066d4d9ed 100644
--- a/dune/gdt/test/pn-discretization.hh
+++ b/dune/gdt/test/pn-discretization.hh
@@ -8,6 +8,8 @@
#ifndef DUNE_GDT_TEST_HYPERBOLIC_PN_DISCRETIZATION_HH
#define DUNE_GDT_TEST_HYPERBOLIC_PN_DISCRETIZATION_HH
+#include <dune/common/exceptions.hh>
+
#include <dune/xt/common/parallel/threadmanager.hh>
#include <dune/xt/common/string.hh>
#include <dune/xt/common/test/gtest/gtest.h>
@@ -33,95 +35,92 @@
#include <dune/gdt/test/momentmodels/kineticequation.hh>
-int parse_momentmodel_arguments(int argc,
- char** argv,
- size_t& num_threads,
- size_t& threading_partition_factor,
- size_t& num_save_steps,
- size_t& num_output_steps,
- size_t& quad_order,
- std::string& grid_size,
- std::string& overlap_size,
- double& t_end,
- std::string& filename)
+void parse_momentmodel_arguments(int argc,
+ char** argv,
+ size_t& num_threads,
+ size_t& threading_partition_factor,
+ size_t& num_save_steps,
+ size_t& num_output_steps,
+ size_t& quad_order,
+ size_t& quad_refinements,
+ size_t& grid_size,
+ size_t& overlap_size,
+ double& t_end,
+ std::string& filename)
{
using namespace Dune;
using namespace Dune::GDT;
MPIHelper::instance(argc, argv);
+ // default values
+ num_threads = 1;
+ threading_partition_factor = 1;
+ num_save_steps = 10;
+ num_output_steps = num_save_steps;
+ quad_order = -1;
+ quad_refinements = -1;
+ grid_size = -1, overlap_size = 2;
+ t_end = 0.;
+ filename = "";
+
for (int i = 1; i < argc; ++i) {
if (std::string(argv[i]) == "--num_threads") {
- if (i + 1 < argc) {
+ if (i + 1 < argc)
num_threads = XT::Common::from_string<size_t>(argv[++i]);
- } else {
- std::cerr << "--num_threads option requires one argument." << std::endl;
- return 1;
- }
+ else
+ DUNE_THROW(Dune::IOError, "--num_threads option requires one argument.");
} else if (std::string(argv[i]) == "--threading_partition_factor") {
- if (i + 1 < argc) {
+ if (i + 1 < argc)
threading_partition_factor = XT::Common::from_string<size_t>(argv[++i]);
- } else {
- std::cerr << "--threading_partition_factor option requires one argument." << std::endl;
- return 1;
- }
+ else
+ DUNE_THROW(Dune::IOError, "--threading_partition_factor option requires one argument.");
} else if (std::string(argv[i]) == "--filename") {
- if (i + 1 < argc) {
+ if (i + 1 < argc)
filename = std::string(argv[++i]);
- } else {
- std::cerr << "--filename option requires one argument." << std::endl;
- return 1;
- }
+ else
+ DUNE_THROW(Dune::IOError, "--filename option requires one argument.");
} else if (std::string(argv[i]) == "--quad_order") {
- if (i + 1 < argc) {
+ if (i + 1 < argc)
quad_order = XT::Common::from_string<size_t>(argv[++i]);
- } else {
- std::cerr << "--quad_order option requires one argument." << std::endl;
- return 1;
- }
+ else
+ DUNE_THROW(Dune::IOError, "--quad_order option requires one argument.");
+ } else if (std::string(argv[i]) == "--quad_refinements") {
+ if (i + 1 < argc)
+ quad_refinements = XT::Common::from_string<size_t>(argv[++i]);
+ else
+ DUNE_THROW(Dune::IOError, "--quad_refinements option requires one argument.");
} else if (std::string(argv[i]) == "--num_save_steps") {
- if (i + 1 < argc) {
+ if (i + 1 < argc)
num_save_steps = XT::Common::from_string<size_t>(argv[++i]);
- } else {
- std::cerr << "--num_save_steps option requires one argument." << std::endl;
- return 1;
- }
+ else
+ DUNE_THROW(Dune::IOError, "--num_save_steps option requires one argument.");
} else if (std::string(argv[i]) == "--num_output_steps") {
- if (i + 1 < argc) {
+ if (i + 1 < argc)
num_output_steps = XT::Common::from_string<size_t>(argv[++i]);
- } else {
- std::cerr << "--num_output_steps option requires one argument." << std::endl;
- return 1;
- }
+ else
+ DUNE_THROW(Dune::IOError, "--num_output_steps option requires one argument.");
} else if (std::string(argv[i]) == "--grid_size") {
- if (i + 1 < argc) {
- grid_size = argv[++i];
- } else {
- std::cerr << "--grid_size option requires one argument." << std::endl;
- return 1;
- }
+ if (i + 1 < argc)
+ grid_size = XT::Common::from_string<size_t>(argv[++i]);
+ else
+ DUNE_THROW(Dune::IOError, "--grid_size option requires one argument.");
} else if (std::string(argv[i]) == "--overlap_size") {
- if (i + 1 < argc) {
- overlap_size = argv[++i];
- } else {
- std::cerr << "--overlap_size option requires one argument." << std::endl;
- return 1;
- }
+ if (i + 1 < argc)
+ overlap_size = XT::Common::from_string<size_t>(argv[++i]);
+ else
+ DUNE_THROW(Dune::IOError, "--overlap_size option requires one argument.");
} else if (std::string(argv[i]) == "--t_end") {
- if (i + 1 < argc) {
+ if (i + 1 < argc)
t_end = XT::Common::from_string<double>(argv[++i]);
- } else {
- std::cerr << "--t_end option requires one argument." << std::endl;
- return 1;
- }
+ else
+ DUNE_THROW(Dune::IOError, "--t_end option requires one argument.");
} else {
- std::cerr << "Unknown option " << std::string(argv[i]) << std::endl;
- return 1;
+ DUNE_THROW(Dune::IOError, "Unknown option " + std::string(argv[i]));
}
}
DXTC_CONFIG.set("threading.partition_factor", threading_partition_factor, true);
XT::Common::threadManager().set_max_threads(num_threads);
- return 0;
}
@@ -176,9 +175,10 @@ struct HyperbolicPnDiscretization
Dune::FieldVector<double, 3> run(size_t num_save_steps = 1,
size_t num_output_steps = 0,
size_t quad_order = size_t(-1),
- std::string grid_size = "",
- std::string overlap_size = "2",
- double t_end = TestCaseType::t_end,
+ size_t quad_refinements = size_t(-1),
+ size_t grid_size = size_t(-1),
+ size_t overlap_size = 2,
+ double t_end = 0.,
std::string filename = "")
{
using namespace Dune;
@@ -205,9 +205,9 @@ struct HyperbolicPnDiscretization
//******************* create grid and FV space ***************************************
auto grid_config = ProblemType::default_grid_cfg();
- if (!grid_size.empty())
- grid_config["num_elements"] = grid_size;
- grid_config["overlap_size"] = overlap_size;
+ if (grid_size != size_t(-1))
+ grid_config["num_elements"] = XT::Common::to_string(grid_size);
+ grid_config["overlap_size"] = XT::Common::to_string(overlap_size);
const auto grid_ptr =
Dune::XT::Grid::CubeGridProviderFactory<GridType>::create(grid_config, MPIHelper::getCommunicator()).grid_ptr();
assert(grid_ptr->comm().size() == 1 || grid_ptr->overlapSize(0) > 0);
@@ -216,9 +216,9 @@ struct HyperbolicPnDiscretization
const AdvectionSourceSpaceType advection_source_space(grid_view, 1);
//******************* create EquationType object ***************************************
- std::shared_ptr<const MomentBasis> basis_functions = (quad_order == size_t(-1))
- ? std::make_shared<const MomentBasis>()
- : std::make_shared<const MomentBasis>(quad_order);
+ std::shared_ptr<const MomentBasis> basis_functions = std::make_shared<const MomentBasis>(
+ quad_order == size_t(-1) ? MomentBasis::default_quad_order() : quad_order,
+ quad_refinements == size_t(-1) ? MomentBasis::default_quad_refinements() : quad_refinements);
const std::unique_ptr<ProblemType> problem_ptr =
XT::Common::make_unique<ProblemType>(*basis_functions, grid_config);
const auto& problem = *problem_ptr;
@@ -307,8 +307,10 @@ struct HyperbolicPnDiscretization
FvOperatorType& fv_operator =
FvOperatorChooser<TestCaseType::reconstruction>::choose(advection_operator, reconstruction_fv_operator);
- filename += "_" + ProblemType::static_id();
- filename += "_grid_" + grid_size;
+ if (!filename.empty())
+ filename += "_";
+ filename += ProblemType::static_id();
+ filename += "_grid_" + grid_config["num_elements"];
filename += "_tend_" + XT::Common::to_string(t_end);
filename += TestCaseType::reconstruction ? "_ord2" : "_ord1";
filename += "_" + basis_functions->pn_name();
@@ -366,7 +368,8 @@ struct HyperbolicPnTest
{
void run()
{
- auto norms = HyperbolicPnDiscretization<TestCaseType>::run();
+ auto norms = HyperbolicPnDiscretization<TestCaseType>::run(
+ 1, 0, size_t(-1), size_t(-1), size_t(-1), 2, TestCaseType::t_end, "test");
const double l1norm = norms[0];
const double l2norm = norms[1];
const double linfnorm = norms[2];