diff --git a/dune/gdt/operators/advection-fv-entropybased.hh b/dune/gdt/operators/advection-fv-entropybased.hh
index 9177813509e4a3a5236adb0895006811700c47dd..44d30970508f0d015b8dcce6b45ba68c6f72dbd8 100644
--- a/dune/gdt/operators/advection-fv-entropybased.hh
+++ b/dune/gdt/operators/advection-fv-entropybased.hh
@@ -11,6 +11,8 @@
#ifndef DUNE_GDT_OPERATORS_FV_ENTROPYBASED_HH
#define DUNE_GDT_OPERATORS_FV_ENTROPYBASED_HH
+#include <boost/config.hpp>
+
#include <dune/xt/common/numeric.hh>
#include <dune/gdt/operators/interfaces.hh>
@@ -158,10 +160,12 @@ public:
static constexpr size_t dimDomain = AdvectionOperatorType::SGV::dimension;
static constexpr size_t dimRange = AdvectionOperatorType::s_r;
static constexpr double interval_length = 2. / (dimRange - 1.);
+ static constexpr size_t first_positive_index_1d = dimRange / 2;
using RangeFieldType = typename BaseType::FieldType;
using ParameterFunctionType = typename ProblemType::ScalarFunctionType;
using DomainType = FieldVector<RangeFieldType, dimDomain>;
using BoundaryFluxesMapType = std::map<DomainType, DynamicVector<RangeFieldType>, XT::Common::FieldVectorFloatLess>;
+ using IntersectionType = XT::Grid::extract_intersection_t<typename AdvectionOperatorType::SGV>;
EntropicCoordinatesMasslumpedOperator(const AdvectionOperatorType& advection_op,
const ProblemType& problem,
@@ -174,6 +178,7 @@ public:
, sigma_s_(problem.sigma_s())
, Q_(problem.Q())
, exp_evaluations_(source_space().mapper().size())
+ , dummy_range_(dimRange)
, dx_(dx)
, h_inv_(1. / dx_)
, quad_points_(dimDomain * dimRange)
@@ -235,48 +240,30 @@ public:
assert(source.size() < std::numeric_limits<int>::max());
XT::Common::Mkl::exp(static_cast<int>(source.size()), source.data(), exp_evaluations_.data());
auto* range_data = range.data();
- std::fill(range_data, range_data + range.size(), 0.);
const auto& grid_view = source_space().grid_view();
const size_t num_entities = grid_view.size(0);
- DomainType entity_center;
if constexpr (dimDomain == 1) {
+ // left boundary flux
+ const auto& left_boundary_flux = boundary_fluxes_.at(lower_left_);
+ for (size_t jj = 0; jj < dimRange; ++jj)
+ range_data[jj] += left_boundary_flux[jj] * h_inv_;
+ // inner fluxes and rhs
for (size_t ii = 0; ii < num_entities; ++ii) {
- const bool left_boundary = (ii == 0);
- const bool right_boundary = (ii == num_entities - 1);
const auto offset = ii * dimRange;
auto* range_entity = range_data + offset;
- auto* range_left = range_entity - dimRange;
- auto* range_right = range_entity + dimRange;
+ auto* range_left = (ii == 0 ? dummy_range_.data() : range_entity - dimRange);
+ auto* range_right = (ii == num_entities - 1 ? dummy_range_.data() : range_entity + dimRange);
const auto* psi = exp_evaluations_.data() + offset;
- const auto* boundary_flux =
- left_boundary ? &(boundary_fluxes_.at(lower_left_)) : &(boundary_fluxes_.at(upper_right_));
- entity_center[0] = lower_left_[0] + (ii + 0.5) * dx_;
- const auto sigma_a_value = sigma_a_->evaluate(entity_center)[0];
- const auto sigma_s_value = sigma_s_->evaluate(entity_center)[0];
- const auto sigma_t_value = sigma_a_value + sigma_s_value;
- const auto Q_value = Q_->evaluate(entity_center)[0];
- auto density = XT::Common::reduce(psi + 1, psi + dimRange - 1, 0.);
- density += 0.5 * (psi[0] + psi[dimRange - 1]);
- density *= interval_length;
- for (size_t jj = 0; jj < dimRange; ++jj) {
- // flux
- const auto& v = quad_points_[jj];
- const auto weight = (jj == 0 || jj == dimRange - 1) ? interval_length / 2. : interval_length;
- const auto flux_jj = std::abs(v) * psi[jj] * h_inv_ * weight;
- range_entity[jj] -= flux_jj;
- if (v > 0.) {
- if (!right_boundary)
- range_right[jj] += flux_jj;
- } else if (!left_boundary) {
- range_left[jj] += flux_jj;
- }
- if (left_boundary || right_boundary)
- range_entity[jj] += (*boundary_flux)[jj] * h_inv_;
- // rhs
- range_entity[jj] +=
- sigma_s_value * density * u_iso_[jj] + Q_value * basis_integrated_[jj] - psi[jj] * weight * sigma_t_value;
- } // jj
+ // flux
+ flux_1d(range_entity, range_left, range_right, psi);
+ // rhs
+ entity_center_[0] = lower_left_[0] + (ii + 0.5) * dx_;
+ rhs_1d(range_entity, psi);
} // entities
+ // left boundary flux
+ const auto& right_boundary_flux = boundary_fluxes_.at(upper_right_);
+ for (size_t jj = 0; jj < dimRange; ++jj)
+ range_data[(num_entities - 1) * dimRange + jj] += right_boundary_flux[jj] * h_inv_;
} else if constexpr (dimDomain == 3) {
for (auto&& entity : Dune::elements(grid_view)) {
const auto entity_index = grid_view.indexSet().index(entity);
@@ -286,54 +273,124 @@ public:
// flux
for (auto&& intersection : Dune::intersections(grid_view, entity)) {
const bool boundary = !intersection.neighbor();
- const auto& outside_entity = boundary ? entity : intersection.outside();
- const auto outside_index = grid_view.indexSet().index(outside_entity);
- auto* range_outside = range_data + outside_index * dimRange;
- const auto dd = intersection.indexInInside() / 2;
- const auto& n = intersection.centerUnitOuterNormal();
- const auto* boundary_flux = boundary ? &(boundary_fluxes_.at(intersection.geometry().center())) : nullptr;
- for (size_t jj = 0; jj < dimRange; ++jj) {
- const auto v_times_n = quad_points_[dd * dimRange + jj] * n[dd];
- if (v_times_n > 0.) {
- const auto flux_jj = v_times_n * psi[jj] * quad_weights_[jj] * h_inv_;
- range_entity[jj] -= flux_jj;
- if (!boundary)
- range_outside[jj] += flux_jj;
- } else if (boundary) {
- range_entity[jj] += (*boundary_flux)[jj] * h_inv_;
- }
- } // jj
+ const auto& outside_entity = intersection.outside();
+ auto* range_outside =
+ boundary ? dummy_range_.data() : range_data + grid_view.indexSet().index(outside_entity) * dimRange;
+ flux_3d(range_entity, range_outside, psi, intersection);
} // intersections
// rhs
- entity_center = entity.geometry().center();
- const auto sigma_a_value = sigma_a_->evaluate(entity_center)[0];
- const auto sigma_s_value = sigma_s_->evaluate(entity_center)[0];
- const auto sigma_t_value = sigma_a_value + sigma_s_value;
- const auto Q_value = Q_->evaluate(entity_center)[0];
- const auto density = XT::Common::transform_reduce(psi, psi + dimRange, quad_weights_.begin(), 0.);
- for (size_t jj = 0; jj < dimRange; ++jj) {
- range_entity[jj] += sigma_s_value * density * u_iso_[jj] + Q_value * basis_integrated_[jj]
- - psi[jj] * quad_weights_[jj] * sigma_t_value;
- } // jj
+ entity_center_ = entity.geometry().center();
+ rhs_3d(range_entity, psi);
} // entities
}
// inverse Hessian
const auto* psi = exp_evaluations_.data();
+ if constexpr (dimDomain == 1)
+ apply_inverse_hessian_1d(num_entities, range_data, psi);
+ else
+ apply_inverse_hessian_3d(num_entities, range_data, psi);
+ } // void apply(...)
+
+ void flux_1d(RangeFieldType* BOOST_RESTRICT range_entity,
+ RangeFieldType* BOOST_RESTRICT range_left,
+ RangeFieldType* BOOST_RESTRICT range_right,
+ const RangeFieldType* BOOST_RESTRICT psi) const
+ {
+ // fluxes in negative direction
+ for (size_t jj = 0; jj < first_positive_index_1d; ++jj) {
+ const auto weight = (jj == 0 ? interval_length * 0.5 : interval_length);
+ const auto flux_jj = -quad_points_[jj] * psi[jj] * h_inv_ * weight;
+ range_entity[jj] -= flux_jj;
+ range_left[jj] += flux_jj;
+ }
+ // fluxes in positive direction
+ for (size_t jj = first_positive_index_1d; jj < dimRange; ++jj) {
+ const auto weight = (jj == dimRange - 1 ? interval_length * 0.5 : interval_length);
+ const auto flux_jj = quad_points_[jj] * psi[jj] * h_inv_ * weight;
+ range_entity[jj] -= flux_jj;
+ range_right[jj] += flux_jj;
+ }
+ }
+
+ void flux_3d(RangeFieldType* BOOST_RESTRICT range_entity,
+ RangeFieldType* BOOST_RESTRICT range_outside,
+ const RangeFieldType* BOOST_RESTRICT psi,
+ const IntersectionType& intersection) const
+ {
+ const bool boundary = !intersection.neighbor();
+ const auto dd = intersection.indexInInside() / 2;
+ const auto n_dd = intersection.centerUnitOuterNormal()[dd];
+ const auto quad_offset = dd * dimRange;
+ for (size_t jj = 0; jj < dimRange; ++jj) {
+ const auto v_times_n = quad_points_[quad_offset + jj] * n_dd;
+ if (v_times_n > 0.) {
+ const auto flux_jj = psi[jj] * quad_weights_[jj] * h_inv_ * v_times_n;
+ range_entity[jj] -= flux_jj;
+ range_outside[jj] += flux_jj;
+ }
+ } // jj
+ if (boundary) {
+ const auto& boundary_flux = boundary_fluxes_.at(intersection.geometry().center());
+ for (size_t jj = 0; jj < dimRange; ++jj)
+ range_entity[jj] += boundary_flux[jj] * h_inv_;
+ }
+ }
+
+ // Note: entity_center_ has to be set to entity.geometry().center() before calling this function
+ void rhs_1d(RangeFieldType* BOOST_RESTRICT range_entity, const RangeFieldType* BOOST_RESTRICT psi) const
+ {
+ const auto sigma_a_value = sigma_a_->evaluate(entity_center_)[0];
+ const auto sigma_s_value = sigma_s_->evaluate(entity_center_)[0];
+ const auto Q_value = Q_->evaluate(entity_center_)[0];
+ auto density = XT::Common::reduce(psi + 1, psi + dimRange - 1, 0.);
+ density += 0.5 * (psi[0] + psi[dimRange - 1]);
+ density *= interval_length;
+ const auto sigma_s_factor = sigma_s_value * density;
+ const auto sigma_t_factor = (sigma_a_value + sigma_s_value) * interval_length;
+ range_entity[0] += u_iso_[0] * sigma_s_factor + basis_integrated_[0] * Q_value - psi[0] * sigma_t_factor * 0.5;
+ for (size_t jj = 1; jj < dimRange - 1; ++jj)
+ range_entity[jj] += u_iso_[jj] * sigma_s_factor + basis_integrated_[jj] * Q_value - psi[jj] * sigma_t_factor;
+ range_entity[dimRange - 1] += u_iso_[dimRange - 1] * sigma_s_factor + basis_integrated_[dimRange - 1] * Q_value
+ - psi[dimRange - 1] * sigma_t_factor * 0.5;
+ }
+
+ // Note: entity_center_ has to be set to entity.geometry().center() before calling this function
+ void rhs_3d(RangeFieldType* BOOST_RESTRICT range_entity, const RangeFieldType* BOOST_RESTRICT psi) const
+ {
+ const auto sigma_a_value = sigma_a_->evaluate(entity_center_)[0];
+ const auto sigma_s_value = sigma_s_->evaluate(entity_center_)[0];
+ const auto sigma_t_value = sigma_a_value + sigma_s_value;
+ const auto Q_value = Q_->evaluate(entity_center_)[0];
+ const auto density = XT::Common::transform_reduce(psi, psi + dimRange, quad_weights_.begin(), 0.);
+ const auto sigma_s_factor = sigma_s_value * density;
+ for (size_t jj = 0; jj < dimRange; ++jj)
+ range_entity[jj] +=
+ u_iso_[jj] * sigma_s_factor + basis_integrated_[jj] * Q_value - psi[jj] * quad_weights_[jj] * sigma_t_value;
+ }
+
+ static void apply_inverse_hessian_1d(const size_t num_entities,
+ RangeFieldType* BOOST_RESTRICT range_data,
+ const RangeFieldType* BOOST_RESTRICT psi)
+ {
for (size_t ii = 0; ii < num_entities; ++ii) {
const auto offset = ii * dimRange;
- for (size_t jj = 0; jj < dimRange; ++jj) {
- auto& val = range_data[offset + jj];
- if constexpr (dimDomain == 1) {
- const auto weight = (jj == 0 || jj == dimRange - 1) ? interval_length / 2. : interval_length;
- val /= psi[offset + jj] * weight;
- } else {
- val /= psi[offset + jj] * quad_weights_[jj];
- }
- if (std::isnan(val) || std::isinf(val))
- DUNE_THROW(Dune::MathError, "inf or nan in range!");
- } // jj
+ range_data[offset] /= psi[offset] * interval_length * 0.5;
+ for (size_t jj = 1; jj < dimRange - 1; ++jj)
+ range_data[offset + jj] /= psi[offset + jj] * interval_length;
+ range_data[offset + dimRange - 1] /= psi[offset + dimRange - 1] * interval_length * 0.5;
} // ii
- } // void apply(...)
+ }
+
+ void apply_inverse_hessian_3d(const size_t num_entities,
+ RangeFieldType* BOOST_RESTRICT range_data,
+ const RangeFieldType* BOOST_RESTRICT psi) const
+ {
+ for (size_t ii = 0; ii < num_entities; ++ii) {
+ const auto offset = ii * dimRange;
+ for (size_t jj = 0; jj < dimRange; ++jj)
+ range_data[offset + jj] /= psi[offset + jj] * quad_weights_[jj];
+ } // ii
+ }
const AdvectionOperatorType& advection_op_;
DynamicVector<RangeFieldType> u_iso_;
@@ -342,6 +399,8 @@ public:
std::unique_ptr<ParameterFunctionType> sigma_s_;
std::unique_ptr<ParameterFunctionType> Q_;
mutable std::vector<RangeFieldType> exp_evaluations_;
+ mutable DomainType entity_center_;
+ mutable std::vector<RangeFieldType> dummy_range_;
const RangeFieldType dx_;
const RangeFieldType h_inv_;
std::vector<RangeFieldType> quad_points_;
@@ -349,7 +408,7 @@ public:
const BoundaryFluxesMapType& boundary_fluxes_;
const DomainType lower_left_;
const DomainType upper_right_;
-}; // class EntropicCoordinatesMasslumpedOperator<...>
+}; // namespace Dune
template <class AdvectionOperatorImp, class EntropySolverImp>
diff --git a/dune/gdt/test/momentmodels/entropyflux_implementations.hh b/dune/gdt/test/momentmodels/entropyflux_implementations.hh
index 012121be623992a16df74fee54f6ca81f9d5cf25..9ef798fc3a1fd41ec57b1a8fb84ce01d964541d5 100644
--- a/dune/gdt/test/momentmodels/entropyflux_implementations.hh
+++ b/dune/gdt/test/momentmodels/entropyflux_implementations.hh
@@ -5277,7 +5277,7 @@ public:
evaluate_eta_ast_prime(alpha_j, eta_ast_prime_vals[1]);
// calculate \sum_{i=1}^d < \omega_i m G_\alpha(u) > n_i
DomainType ret(0);
- const size_t first_positive = dimRange / 2;
+ constexpr size_t first_positive = dimRange / 2;
const auto& first_exp_vals = n_ij[dd] < 0. ? eta_ast_prime_vals[0] : eta_ast_prime_vals[1];
for (size_t ll = 0; ll < first_positive; ++ll)
ret[ll] = first_exp_vals[ll] * grid_points_[ll];
@@ -5296,7 +5296,7 @@ public:
evaluate_eta_ast_prime(alpha_i, eta_ast_prime_vals);
// calculate \sum_{i=1}^d < \omega_i m G_\alpha(u) > n_i
DomainType ret(0);
- const size_t first_positive = dimRange / 2;
+ constexpr size_t first_positive = dimRange / 2;
if (n_ij[dd] < 0.) {
for (size_t ll = 0; ll < first_positive; ++ll)
ret[ll] = eta_ast_prime_vals[ll] * grid_points_[ll];
@@ -5334,7 +5334,7 @@ public:
const auto& psi_right = (*ansatz_distribution_values[2]);
constexpr bool reconstruct = (slope_type != SlopeLimiterType::no_slope);
// reconstruct densities
- const size_t first_positive = dimRange / 2;
+ constexpr size_t first_positive = dimRange / 2;
for (size_t ll = 0; ll < first_positive; ++ll) {
const RangeFieldType slope =
reconstruct ? slope_func(psi_entity[ll] - psi_left[ll], psi_right[ll] - psi_entity[ll]) : 0.;
diff --git a/dune/gdt/test/momentmodels/kinetictransport/checkerboard.hh b/dune/gdt/test/momentmodels/kinetictransport/checkerboard.hh
index 4f90aa9a75323873558ac3021d083cc9469435a5..9e9916a1126d5c4dd8fb995ee891f463c43e8546 100644
--- a/dune/gdt/test/momentmodels/kinetictransport/checkerboard.hh
+++ b/dune/gdt/test/momentmodels/kinetictransport/checkerboard.hh
@@ -102,10 +102,10 @@ public:
protected:
static bool is_absorbing(const FieldVector<RangeFieldType, 2>& x)
{
- size_t row = XT::Common::numeric_cast<size_t>(std::floor(x[1]));
+ size_t row = static_cast<size_t>(std::floor(x[1]));
if (row == 7)
row = 6;
- size_t col = XT::Common::numeric_cast<size_t>(std::floor(x[0]));
+ size_t col = static_cast<size_t>(std::floor(x[0]));
if (col == 7)
col = 6;
assert(row < 7 && col < 7);
@@ -115,9 +115,9 @@ protected:
static bool is_absorbing(const FieldVector<RangeFieldType, 3>& x)
{
- size_t plane = XT::Common::numeric_cast<size_t>(std::floor(x[2]));
- size_t col = XT::Common::numeric_cast<size_t>(std::floor(x[0]));
- size_t row = XT::Common::numeric_cast<size_t>(std::floor(x[1]));
+ size_t plane = static_cast<size_t>(std::floor(x[2]));
+ size_t col = static_cast<size_t>(std::floor(x[0]));
+ size_t row = static_cast<size_t>(std::floor(x[1]));
assert(plane <= 7 && row <= 7 && col <= 7);
if (plane == 0 || plane >= 6)
return false;
@@ -130,16 +130,16 @@ protected:
static bool is_center(const FieldVector<RangeFieldType, 2>& x)
{
- size_t row = XT::Common::numeric_cast<size_t>(std::floor(x[1]));
- size_t col = XT::Common::numeric_cast<size_t>(std::floor(x[0]));
+ size_t row = static_cast<size_t>(std::floor(x[1]));
+ size_t col = static_cast<size_t>(std::floor(x[0]));
return row == 3 && col == 3;
}
static bool is_center(const FieldVector<RangeFieldType, 3>& x)
{
- size_t plane = XT::Common::numeric_cast<size_t>(std::floor(x[2]));
- size_t row = XT::Common::numeric_cast<size_t>(std::floor(x[1]));
- size_t col = XT::Common::numeric_cast<size_t>(std::floor(x[0]));
+ size_t plane = static_cast<size_t>(std::floor(x[2]));
+ size_t row = static_cast<size_t>(std::floor(x[1]));
+ size_t col = static_cast<size_t>(std::floor(x[0]));
return plane == 3 && row == 3 && col == 3;
}
}; // class CheckerboardPn<...>
diff --git a/dune/gdt/tools/timestepper/adaptive-rungekutta-kinetic.hh b/dune/gdt/tools/timestepper/adaptive-rungekutta-kinetic.hh
index b8806d268cb797cd1f2afe47646bea564d8f4e6d..b3e0be8cb6bdb93985d31b46466c90dc39378904 100644
--- a/dune/gdt/tools/timestepper/adaptive-rungekutta-kinetic.hh
+++ b/dune/gdt/tools/timestepper/adaptive-rungekutta-kinetic.hh
@@ -200,20 +200,22 @@ public:
auto& t = current_time();
auto& alpha_n = current_solution();
+ const auto num_dofs = alpha_n.dofs().vector().size();
size_t first_stage_to_compute = 0;
if (first_same_as_last_ && last_stage_of_previous_step_) {
stages_k_[0].dofs().vector() = last_stage_of_previous_step_->dofs().vector();
first_stage_to_compute = 1;
}
first_same_as_last_ = true;
- while (mixed_error > 1.) {
+ while (!(mixed_error < 1.)) {
bool skip_error_computation = false;
actual_dt *= time_step_scale_factor;
for (size_t ii = first_stage_to_compute; ii < num_stages_ - 1; ++ii) {
- set_op_param("t", t + actual_dt * c_[ii]), stages_k_[ii].dofs().vector() *= 0.;
+ set_op_param("t", t + actual_dt * c_[ii]);
+ std::fill_n(&(stages_k_[ii].dofs().vector()[0]), num_dofs, 0.);
alpha_tmp_.dofs().vector() = alpha_n.dofs().vector();
for (size_t jj = 0; jj < ii; ++jj)
- alpha_tmp_.dofs().vector().axpy(actual_dt * r_ * (A_[ii][jj]), stages_k_[jj].dofs().vector());
+ alpha_tmp_.dofs().vector().axpy(actual_dt * r_ * A_[ii][jj], stages_k_[jj].dofs().vector());
try {
op_.apply(alpha_tmp_.dofs().vector(), stages_k_[ii].dofs().vector(), op_param_);
if (regularize_if_needed(consider_regularization, r_it, r_sequence)) {
@@ -247,7 +249,7 @@ public:
// calculate last stage
set_op_param("t", t + actual_dt * c_[num_stages_ - 1]);
- stages_k_[num_stages_ - 1].dofs().vector() *= 0.;
+ std::fill_n(&(stages_k_[num_stages_ - 1].dofs().vector()[0]), num_dofs, 0.);
try {
op_.apply(alpha_np1_.dofs().vector(), stages_k_[num_stages_ - 1].dofs().vector(), op_param_);
if (regularize_if_needed(consider_regularization, r_it, r_sequence)) {
@@ -273,31 +275,40 @@ public:
for (size_t ii = 0; ii < num_stages_; ++ii)
alpha_tmp_.dofs().vector().axpy(actual_dt * r_ * b_2_[ii], stages_k_[ii].dofs().vector());
- // calculate error
const auto* alpha_tmp_data =
XT::Common::VectorAbstraction<typename DiscreteFunctionType::VectorType>::data(alpha_tmp_.dofs().vector());
const auto* alpha_np1_data =
XT::Common::VectorAbstraction<typename DiscreteFunctionType::VectorType>::data(alpha_np1_.dofs().vector());
- mixed_error =
- XT::Common::transform_reduce(alpha_tmp_data,
- alpha_tmp_data + alpha_tmp_.dofs().vector().size(),
- alpha_np1_data,
- 0.,
- /*reduction*/ [](const auto& a, const auto& b) { return std::max(a, b); },
- /*transformation*/
- [atol = atol_, rtol = rtol_](const auto& a, const auto& b) {
- return std::abs(a - b) / (atol + std::max(std::abs(a), std::abs(b)) * rtol);
- });
- // std::cout << mixed_error << std::endl;
- // scale dt to get the estimated optimal time step length
- time_step_scale_factor =
- std::min(std::max(0.8 * std::pow(1. / mixed_error, 1. / (q + 1.)), scale_factor_min_), scale_factor_max_);
-
- // maybe adjust alpha to enforce a minimum density or avoid problems with matrix conditions
- if (!(mixed_error > 1.)
- && min_density_setter_.apply_with_dt(alpha_np1_.dofs().vector(), alpha_np1_.dofs().vector(), actual_dt)) {
- // we cannot use the first-same-as-last property for the next time step if we changed alpha
- first_same_as_last_ = false;
+
+ // if b == NaN, std::max(a, b) might return a, so mixed_error might be non-NaN though there are NaNs in the
+ // vectors So check for NaNs before
+ bool nan_found = check_for_nan(alpha_tmp_data, alpha_np1_data, num_dofs);
+ if (nan_found) {
+ mixed_error = 1e10;
+ time_step_scale_factor = 0.5;
+ } else {
+ // calculate error
+ mixed_error = XT::Common::transform_reduce(
+ alpha_tmp_data,
+ alpha_tmp_data + num_dofs,
+ alpha_np1_data,
+ 0.,
+ /*reduction*/ [](const auto& a, const auto& b) { return std::max(a, b); },
+ /*transformation*/
+ [atol = atol_, rtol = rtol_](const auto& a, const auto& b) {
+ return std::abs(a - b) / (atol + std::max(std::abs(a), std::abs(b)) * rtol);
+ });
+ // std::cout << mixed_error << std::endl;
+ // scale dt to get the estimated optimal time step length
+ time_step_scale_factor =
+ std::min(std::max(0.8 * std::pow(1. / mixed_error, 1. / (q + 1.)), scale_factor_min_), scale_factor_max_);
+
+ // maybe adjust alpha to enforce a minimum density or avoid problems with matrix conditions
+ if (mixed_error < 1.
+ && min_density_setter_.apply_with_dt(alpha_np1_.dofs().vector(), alpha_np1_.dofs().vector(), actual_dt)) {
+ // we cannot use the first-same-as-last property for the next time step if we changed alpha
+ first_same_as_last_ = false;
+ }
}
} // if (!skip_error_computation)
} // while (mixed_error > 1.)
@@ -315,6 +326,14 @@ public:
} // ... step(...)
private:
+ bool check_for_nan(const RangeFieldType* vec1, const RangeFieldType* vec2, const size_t num_dofs)
+ {
+ for (size_t ii = 0; ii < num_dofs; ++ii)
+ if (std::isnan(vec1[ii] + vec2[ii]))
+ return true;
+ return false;
+ }
+
const MinDensitySetterType min_density_setter_;
const OperatorType& op_;
const EntropyFluxType& entropy_flux_;