From b2f929202cc3302f7026cf041fe0ebcc133f0d7a Mon Sep 17 00:00:00 2001
From: Tobias Leibner <tobias.leibner@googlemail.com>
Date: Fri, 20 Dec 2019 17:33:13 +0100
Subject: [PATCH] [local.operators.advection-fv] add dynamic evaluate
---
dune/gdt/local/operators/advection-fv.hh | 124 ++++++++++++------
.../test/inviscid-compressible-flow/base.hh | 19 +--
.../entropic-coords-mn-discretization.hh | 22 ++--
dune/gdt/test/momentmodels/entropyflux.hh | 17 ++-
.../momentmodels/entropyflux_kineticcoords.hh | 69 +++++++---
.../test/momentmodels/mn-discretization.hh | 15 ++-
.../test/momentmodels/pn-discretization.hh | 14 +-
7 files changed, 183 insertions(+), 97 deletions(-)
diff --git a/dune/gdt/local/operators/advection-fv.hh b/dune/gdt/local/operators/advection-fv.hh
index 51ff9b74b..fdebe9549 100644
--- a/dune/gdt/local/operators/advection-fv.hh
+++ b/dune/gdt/local/operators/advection-fv.hh
@@ -60,7 +60,7 @@ public:
using typename BaseType::LocalSourceType;
using typename BaseType::SourceSpaceType;
using typename BaseType::SourceType;
- using StateType = typename XT::Functions::RangeTypeSelector<SR, m, 1>::type;
+ using DynamicStateType = typename XT::Functions::RangeTypeSelector<SR, m, 1>::dynamic_type;
using NumericalFluxType = NumericalFluxInterface<I, d, m, RR>;
using LocalIntersectionCoords = typename NumericalFluxType::LocalIntersectionCoords;
static const typename LocalSourceType::DomainType static_x;
@@ -71,6 +71,9 @@ public:
: BaseType(2, numerical_flux.parameter_type())
, numerical_flux_(numerical_flux.copy())
, source_is_elementwise_constant_(source_is_elementwise_constant)
+ , u_(m)
+ , v_(m)
+ , g_(m)
{}
LocalAdvectionFvCouplingOperator(const SourceType& source,
@@ -79,6 +82,9 @@ public:
: BaseType(source, 2, numerical_flux.parameter_type())
, numerical_flux_(numerical_flux.copy())
, source_is_elementwise_constant_(source_is_elementwise_constant)
+ , u_(m)
+ , v_(m)
+ , g_(m)
{}
LocalAdvectionFvCouplingOperator(const SourceSpaceType& source_space,
@@ -88,6 +94,9 @@ public:
: BaseType(source_space, source_vector, 2, numerical_flux.parameter_type())
, numerical_flux_(numerical_flux.copy())
, source_is_elementwise_constant_(source_is_elementwise_constant)
+ , u_(m)
+ , v_(m)
+ , g_(m)
{}
LocalAdvectionFvCouplingOperator(const DiscreteSourceType& source, const NumericalFluxType& numerical_flux)
@@ -98,6 +107,9 @@ public:
: BaseType(other)
, numerical_flux_(other.numerical_flux_->copy())
, source_is_elementwise_constant_(other.source_is_elementwise_constant_)
+ , u_(m)
+ , v_(m)
+ , g_(m)
{}
std::unique_ptr<BaseType> copy() const override final
@@ -120,20 +132,23 @@ public:
|| (local_range_outside.space().type() != SpaceType::finite_volume),
Exceptions::operator_error,
"Use LocalAdvectionDgCouplingOperator instead!");
- u_ = local_sources_[0]->evaluate(source_is_elementwise_constant_ ? static_x
- : intersection().geometryInInside().center());
- v_ = local_sources_[1]->evaluate(source_is_elementwise_constant_ ? static_x
- : intersection().geometryInOutside().center());
+ local_sources_[0]->evaluate(
+ source_is_elementwise_constant_ ? static_x : intersection().geometryInInside().center(), u_, param);
+ local_sources_[1]->evaluate(
+ source_is_elementwise_constant_ ? static_x : intersection().geometryInOutside().center(), v_, param);
const auto normal = intersection().centerUnitOuterNormal();
if (numerical_flux_->x_dependent())
x_in_intersection_coords_ = intersection().geometry().local(intersection().geometry().center());
- const auto g = numerical_flux_->apply(x_in_intersection_coords_, u_, v_, normal, param);
+ numerical_flux_->apply(x_in_intersection_coords_, u_, v_, normal, g_, param);
const auto h_intersection = intersection().geometry().volume();
- const auto h_inside_element = intersection().inside().geometry().volume();
- const auto h_outside_element = intersection().outside().geometry().volume();
+ const auto hinv_inside_element = 1. / intersection().inside().geometry().volume();
+ const auto hinv_outside_element = 1. / intersection().outside().geometry().volume();
+ auto& local_range_inside_dofs = local_range_inside.dofs();
+ auto& local_range_outside_dofs = local_range_outside.dofs();
for (size_t ii = 0; ii < m; ++ii) {
- local_range_inside.dofs()[ii] += (g[ii] * h_intersection) / h_inside_element;
- local_range_outside.dofs()[ii] -= (g[ii] * h_intersection) / h_outside_element;
+ const auto g_ii = g_[ii] * h_intersection;
+ local_range_inside_dofs[ii] += g_ii * hinv_inside_element;
+ local_range_outside_dofs[ii] -= g_ii * hinv_outside_element;
}
} // ... apply(...)
@@ -149,8 +164,9 @@ private:
std::unique_ptr<NumericalFluxType> numerical_flux_;
const bool source_is_elementwise_constant_;
mutable LocalIntersectionCoords x_in_intersection_coords_;
- mutable StateType u_;
- mutable StateType v_;
+ mutable DynamicStateType u_;
+ mutable DynamicStateType v_;
+ mutable DynamicStateType g_;
}; // class LocalAdvectionFvCouplingOperator
template <class I,
@@ -186,9 +202,11 @@ public:
using typename BaseType::SourceType;
using StateDomainType = FieldVector<typename SGV::ctype, SGV::dimension>;
- using StateType = typename CouplingOperatorType::StateType;
- using LambdaType = std::function<StateType(
- const StateType& /*u*/, const StateDomainType& /*n*/, const XT::Common::Parameter& /*param*/)>;
+ using DynamicStateType = typename CouplingOperatorType::DynamicStateType;
+ using LambdaType = std::function<void(const DynamicStateType& /*u*/,
+ const StateDomainType& /*n*/,
+ DynamicStateType& /*g*/,
+ const XT::Common::Parameter& /*param*/)>;
// When using this constructor, source has to be set by a call to with_source before calling apply
LocalAdvectionFvBoundaryTreatmentByCustomNumericalFluxOperator(
@@ -198,6 +216,8 @@ public:
: BaseType(1, boundary_treatment_param_type)
, numerical_boundary_flux_(numerical_boundary_flux_lambda)
, source_is_elementwise_constant_(source_is_elementwise_constant)
+ , u_(m)
+ , g_(m)
{}
LocalAdvectionFvBoundaryTreatmentByCustomNumericalFluxOperator(
@@ -208,6 +228,8 @@ public:
: BaseType(source, 1, boundary_treatment_param_type)
, numerical_boundary_flux_(numerical_boundary_flux_lambda)
, source_is_elementwise_constant_(source_is_elementwise_constant)
+ , u_(m)
+ , g_(m)
{}
LocalAdvectionFvBoundaryTreatmentByCustomNumericalFluxOperator(
@@ -219,6 +241,8 @@ public:
: BaseType(source_space, source_vector, 1, boundary_treatment_param_type)
, numerical_boundary_flux_(numerical_boundary_flux_lambda)
, source_is_elementwise_constant_(source_is_elementwise_constant)
+ , u_(m)
+ , g_(m)
{}
LocalAdvectionFvBoundaryTreatmentByCustomNumericalFluxOperator(
@@ -235,6 +259,8 @@ public:
: BaseType(other)
, numerical_boundary_flux_(other.numerical_boundary_flux_)
, source_is_elementwise_constant_(other.source_is_elementwise_constant_)
+ , u_(m)
+ , g_(m)
{}
std::unique_ptr<BaseType> copy() const override final
@@ -257,21 +283,23 @@ public:
DUNE_THROW_IF(local_range_inside.space().type() != SpaceType::finite_volume,
Exceptions::operator_error,
"Use LocalAdvectionDgBoundaryTreatmentByCustomNumericalFluxOperator instead!");
- u_ = local_sources_[0]->evaluate(source_is_elementwise_constant_ ? CouplingOperatorType::static_x
- : intersection().geometryInInside().center());
+ local_sources_[0]->evaluate(source_is_elementwise_constant_ ? CouplingOperatorType::static_x
+ : intersection().geometryInInside().center(),
+ u_);
const auto normal = intersection().centerUnitOuterNormal();
- const auto g = numerical_boundary_flux_(u_, normal, param);
- const auto h_intersection = intersection().geometry().volume();
- const auto h_element = intersection().inside().geometry().volume();
+ numerical_boundary_flux_(u_, normal, g_, param);
+ auto& local_range_inside_dofs = local_range_inside.dofs();
+ const auto factor = intersection().geometry().volume() / intersection().inside().geometry().volume();
for (size_t ii = 0; ii < m; ++ii)
- local_range_inside.dofs()[ii] += (g[ii] * h_intersection) / h_element;
+ local_range_inside_dofs[ii] += g_[ii] * factor;
} // ... apply(...)
private:
using BaseType::local_sources_;
const LambdaType numerical_boundary_flux_;
const bool source_is_elementwise_constant_;
- mutable StateType u_;
+ mutable DynamicStateType u_;
+ mutable DynamicStateType g_;
}; // class LocalAdvectionFvBoundaryTreatmentByCustomNumericalFluxOperator
@@ -296,12 +324,13 @@ public:
using NumericalFluxType = NumericalFluxInterface<I, d, m, RF>;
using LocalIntersectionCoords = typename NumericalFluxType::LocalIntersectionCoords;
using FluxType = typename NumericalFluxType::FluxType;
- using StateType = typename CouplingOperatorType::StateType;
- using LambdaType = std::function<StateType(const IntersectionType& /*intersection*/,
- const FieldVector<D, d - 1>& /*xx_in_reference_intersection_coordinates*/,
- const FluxType& /*flux*/,
- const StateType& /*u*/,
- const XT::Common::Parameter& /*param*/)>;
+ using DynamicStateType = typename CouplingOperatorType::DynamicStateType;
+ using LambdaType = std::function<void(const IntersectionType& /*intersection*/,
+ const FieldVector<D, d - 1>& /*xx_in_reference_intersection_coordinates*/,
+ const FluxType& /*flux*/,
+ const DynamicStateType& /*u*/,
+ DynamicStateType& /*v*/,
+ const XT::Common::Parameter& /*param*/)>;
// When using this constructor, source has to be set by a call to with_source before calling apply
LocalAdvectionFvBoundaryTreatmentByCustomExtrapolationOperator(
@@ -313,6 +342,9 @@ public:
, numerical_flux_(numerical_flux.copy())
, extrapolate_(boundary_extrapolation_lambda)
, source_is_elementwise_constant_(source_is_elementwise_constant)
+ , u_(m)
+ , v_(m)
+ , g_(m)
{}
LocalAdvectionFvBoundaryTreatmentByCustomExtrapolationOperator(
@@ -325,6 +357,9 @@ public:
, numerical_flux_(numerical_flux.copy())
, extrapolate_(boundary_extrapolation_lambda)
, source_is_elementwise_constant_(source_is_elementwise_constant)
+ , u_(m)
+ , v_(m)
+ , g_(m)
{}
LocalAdvectionFvBoundaryTreatmentByCustomExtrapolationOperator(
@@ -357,6 +392,9 @@ public:
, numerical_flux_(other.numerical_flux_->copy())
, extrapolate_(other.extrapolate_)
, source_is_elementwise_constant_(other.source_is_elementwise_constant_)
+ , u_(m)
+ , v_(m)
+ , g_(m)
{}
std::unique_ptr<BaseType> copy() const override final
@@ -380,19 +418,22 @@ public:
"Use LocalAdvectionDgBoundaryTreatmentByCustomExtrapolationOperator instead!");
if (numerical_flux_->x_dependent())
x_in_intersection_coords_ = intersection().geometry().local(intersection().geometry().center());
- u_ = local_sources_[0]->evaluate(source_is_elementwise_constant_ ? CouplingOperatorType::static_x
- : intersection().geometryInInside().center());
- v_ = extrapolate_(intersection(),
- ReferenceElements<D, d - 1>::general(intersection().type()).position(0, 0),
- numerical_flux_->flux(),
- u_,
- param);
+ local_sources_[0]->evaluate(source_is_elementwise_constant_ ? CouplingOperatorType::static_x
+ : intersection().geometryInInside().center(),
+ u_,
+ param);
+ extrapolate_(intersection(),
+ ReferenceElements<D, d - 1>::general(intersection().type()).position(0, 0),
+ numerical_flux_->flux(),
+ u_,
+ v_,
+ param);
const auto normal = intersection().centerUnitOuterNormal();
- const auto g = numerical_flux_->apply(x_in_intersection_coords_, u_, v_, normal, param);
- const auto h_intersection = intersection().geometry().volume();
- const auto h_element = intersection().inside().geometry().volume();
+ numerical_flux_->apply(x_in_intersection_coords_, u_, v_, normal, g_, param);
+ auto& local_range_inside_dofs = local_range_inside.dofs();
+ const auto factor = intersection().geometry().volume() / intersection().inside().geometry().volume();
for (size_t ii = 0; ii < m; ++ii)
- local_range_inside.dofs()[ii] += (g[ii] * h_intersection) / h_element;
+ local_range_inside_dofs[ii] += g_[ii] * factor;
} // ... apply(...)
protected:
@@ -408,8 +449,9 @@ private:
const LambdaType extrapolate_;
const bool source_is_elementwise_constant_;
mutable LocalIntersectionCoords x_in_intersection_coords_;
- mutable StateType u_;
- mutable StateType v_;
+ mutable DynamicStateType u_;
+ mutable DynamicStateType v_;
+ mutable DynamicStateType g_;
}; // class LocalAdvectionFvBoundaryTreatmentByCustomExtrapolationOperator
diff --git a/dune/gdt/test/inviscid-compressible-flow/base.hh b/dune/gdt/test/inviscid-compressible-flow/base.hh
index 2f8a45915..31f63ade3 100644
--- a/dune/gdt/test/inviscid-compressible-flow/base.hh
+++ b/dune/gdt/test/inviscid-compressible-flow/base.hh
@@ -195,8 +195,8 @@ protected:
/*periodicity_restriction=*/impermeable_wall_filter);
// the actual handling of impermeable walls
op->append(/*numerical_boundary_flux=*/
- [&](const auto& u, const auto& n, const auto& /*param*/) {
- return euler_tools.flux_at_impermeable_walls(XT::LA::convert_to<RangeType>(u), n);
+ [&](const auto& u, const auto& n, auto& ret, const auto& /*param*/) {
+ ret = euler_tools.flux_at_impermeable_walls(XT::LA::convert_to<RangeType>(u), n);
},
{},
impermeable_wall_filter);
@@ -219,6 +219,7 @@ protected:
const auto& xx_in_reference_intersection_coordinates,
const auto& /*flux*/,
const auto& u,
+ auto& v,
const auto& /*param*/) {
const auto normal = intersection.unitOuterNormal(xx_in_reference_intersection_coordinates);
const auto rho_v_p = euler_tools.primitives(XT::LA::convert_to<RangeType>(u));
@@ -226,7 +227,7 @@ protected:
auto velocity = std::get<1>(rho_v_p);
const auto& pressure = std::get<2>(rho_v_p);
velocity -= normal * 2. * (velocity * normal);
- return euler_tools.conservative(rho, velocity, pressure);
+ v = euler_tools.conservative(rho, velocity, pressure);
},
{},
impermeable_wall_filter);
@@ -264,6 +265,7 @@ protected:
const auto& xx_in_reference_intersection_coordinates,
const auto& /*flux*/,
const auto& u_,
+ auto& v,
const auto& param) {
// evaluate boundary values
const auto element = intersection.inside();
@@ -280,22 +282,22 @@ protected:
// compute v
if (flow_speed < -a) {
// supersonic inlet
- return bv;
+ v = bv;
} else if (!(flow_speed > 0)) {
// subsonic inlet
const auto rho_outer = euler_tools.density(bv);
const auto v_outer = euler_tools.velocity(bv);
const auto p_inner = euler_tools.pressure(u);
- return euler_tools.conservative(rho_outer, v_outer, p_inner);
+ v = euler_tools.conservative(rho_outer, v_outer, p_inner);
} else if (flow_speed < a) {
// subsonic outlet
const auto rho_inner = euler_tools.density(u);
const auto v_inner = euler_tools.velocity(u);
const auto p_outer = euler_tools.pressure(bv);
- return euler_tools.conservative(rho_inner, v_inner, p_outer);
+ v = euler_tools.conservative(rho_inner, v_inner, p_outer);
} else {
// supersonic outlet
- return RangeType(u);
+ v = u;
}
}; // ... heuristic_euler_inflow_outflow_treatment(...)
op->append(heuristic_euler_inflow_outflow_treatment, {}, inflow_outflow_filter);
@@ -306,6 +308,7 @@ protected:
const auto& xx_in_reference_intersection_coordinates,
const auto& /*flux*/,
const auto& u,
+ auto& v,
const auto& /*param*/) {
const auto normal = intersection.unitOuterNormal(xx_in_reference_intersection_coordinates);
const auto rho_v_p = euler_tools.primitives(XT::LA::convert_to<RangeType>(u));
@@ -313,7 +316,7 @@ protected:
auto velocity = std::get<1>(rho_v_p);
const auto& pressure = std::get<2>(rho_v_p);
velocity -= normal * 2. * (velocity * normal);
- return euler_tools.conservative(rho, velocity, pressure);
+ v = euler_tools.conservative(rho, velocity, pressure);
},
{},
impermeable_wall_filter);
diff --git a/dune/gdt/test/momentmodels/entropic-coords-mn-discretization.hh b/dune/gdt/test/momentmodels/entropic-coords-mn-discretization.hh
index cbf2cff91..e547c7619 100644
--- a/dune/gdt/test/momentmodels/entropic-coords-mn-discretization.hh
+++ b/dune/gdt/test/momentmodels/entropic-coords-mn-discretization.hh
@@ -79,6 +79,7 @@ struct HyperbolicEntropicCoordsMnDiscretization
using GenericFunctionType = XT::Functions::GenericFunction<dimDomain, dimRange, 1, RangeFieldType>;
using DomainType = FieldVector<RangeFieldType, dimDomain>;
using RangeType = FieldVector<RangeFieldType, dimRange>;
+ using DynamicRangeType = DynamicVector<RangeFieldType>;
// static const RangeFieldType scale_factor = 1e2;
//******************* create grid and FV space ***************************************
@@ -124,7 +125,9 @@ struct HyperbolicEntropicCoordsMnDiscretization
alpha_boundary_vals.insert(std::make_pair(x, analytical_flux->get_alpha(u)));
}
GenericFunctionType boundary_values_alpha(
- 1, [&](const DomainType& x, const XT::Common::Parameter&) { return alpha_boundary_vals[x]; });
+ 1, [&](const DomainType& x, DynamicRangeType& ret, const XT::Common::Parameter&) {
+ ret = alpha_boundary_vals[x];
+ });
// ***************** project initial values to discrete function *********************
@@ -224,7 +227,6 @@ struct HyperbolicEntropicCoordsMnDiscretization
using BoundaryOperator =
LocalAdvectionFvBoundaryTreatmentByCustomExtrapolationOperator<I, VectorType, GV, dimRange>;
using LambdaType = typename BoundaryOperator::LambdaType;
- using StateType = typename BoundaryOperator::StateType;
// calculate boundary kinetic fluxes
// apply density_operator first to store boundary_evaluations
@@ -233,27 +235,29 @@ struct HyperbolicEntropicCoordsMnDiscretization
density_operator.apply(alpha.dofs().vector(), alpha.dofs().vector());
// store boundary fluxes
- std::map<DomainType, RangeType, XT::Common::FieldVectorFloatLess> boundary_fluxes;
+ std::map<DomainType, DynamicRangeType, 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);
+ const DynamicRangeType 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]; });
+ 1, [&](const DomainType& x, DynamicRangeType& ret, const XT::Common::Parameter&) { ret = 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 boundary_kinetic_fluxes.evaluate(
- intersection.geometry().global(xx_in_reference_intersection_coordinates));
+ const DynamicRangeType& /*u*/,
+ DynamicRangeType& v,
+ const XT::Common::Parameter& param) {
+ boundary_kinetic_fluxes.evaluate(
+ intersection.geometry().global(xx_in_reference_intersection_coordinates), v, param);
};
XT::Grid::ApplyOn::NonPeriodicBoundaryIntersections<GV> filter;
advection_operator.append(boundary_lambda, {}, filter);
diff --git a/dune/gdt/test/momentmodels/entropyflux.hh b/dune/gdt/test/momentmodels/entropyflux.hh
index e10a013de..bbd3e6f7e 100644
--- a/dune/gdt/test/momentmodels/entropyflux.hh
+++ b/dune/gdt/test/momentmodels/entropyflux.hh
@@ -313,12 +313,14 @@ public:
index_set_, entity_caches_, use_thread_cache_, use_entity_cache_, mutexes_, *implementation_);
}
- StateType evaluate_kinetic_flux(const E& inside_entity,
- const E& outside_entity,
- const StateType& u_i,
- const StateType& u_j,
- const DomainType& n_ij,
- const size_t dd) const
+ template <class StateTp, class RetType>
+ void evaluate_kinetic_flux(const E& inside_entity,
+ const E& outside_entity,
+ const StateTp& u_i,
+ const StateTp& u_j,
+ const DomainType& n_ij,
+ const size_t dd,
+ RetType& ret) const
{
// calculate \sum_{i=1}^d < \omega_i m G_\alpha(u) > n_i
const auto local_func = derived_local_function();
@@ -326,9 +328,10 @@ public:
const auto alpha_i = local_func->get_alpha(u_i, true)->first;
local_func->bind(outside_entity);
const auto alpha_j = local_func->get_alpha(u_j, true)->first;
- return implementation_->evaluate_kinetic_flux_with_alphas(alpha_i, alpha_j, n_ij, dd);
+ ret = implementation_->evaluate_kinetic_flux_with_alphas(alpha_i, alpha_j, n_ij, dd);
} // StateType evaluate_kinetic_flux(...)
+
// Returns alpha(u), starting from alpha_iso. To get better performance when calculating several alphas, use
// Localfunction's get_alpha
std::unique_ptr<AlphaReturnType> get_alpha(const StateType& u, const bool regularize) const
diff --git a/dune/gdt/test/momentmodels/entropyflux_kineticcoords.hh b/dune/gdt/test/momentmodels/entropyflux_kineticcoords.hh
index a4ecb753c..8d1b1462c 100644
--- a/dune/gdt/test/momentmodels/entropyflux_kineticcoords.hh
+++ b/dune/gdt/test/momentmodels/entropyflux_kineticcoords.hh
@@ -14,6 +14,7 @@
#include <memory>
#include <dune/xt/common/float_cmp.hh>
+#include <dune/xt/common/numeric.hh>
#include <dune/xt/common/vector_less.hh>
#include <dune/xt/functions/interfaces/flux-function.hh>
@@ -48,6 +49,7 @@ public:
static const size_t dimFlux = MomentBasis::dimFlux;
static const size_t basis_dimRange = MomentBasis::dimRange;
using typename BaseType::DomainType;
+ using typename BaseType::DynamicStateType;
using typename BaseType::E;
using typename BaseType::LocalFunctionType;
using typename BaseType::RangeFieldType;
@@ -119,6 +121,8 @@ public:
const ImplementationType& implementation_;
}; // class Localfunction
+ using RangeReturnType = typename Localfunction::RangeReturnType;
+
bool x_dependent() const override final
{
return false;
@@ -134,14 +138,18 @@ public:
return std::make_unique<Localfunction>(*implementation_);
}
- StateType evaluate_kinetic_flux(const E& /*inside_entity*/,
- const E& /*outside_entity*/,
- const StateType& flux_i,
- const StateType& flux_j,
- const DomainType& n_ij,
- const size_t dd) const
+ template <class StateTp, class RetType>
+ void evaluate_kinetic_flux(const E& /*inside_entity*/,
+ const E& /*outside_entity*/,
+ const StateTp& flux_i,
+ const StateTp& flux_j,
+ const DomainType& n_ij,
+ const size_t dd,
+ RetType& ret) const
{
- return (flux_i + flux_j) * n_ij[dd];
+ ret = flux_i;
+ ret += flux_j;
+ ret *= n_ij[dd];
} // StateType evaluate_kinetic_flux(...)
const MomentBasis& basis_functions() const
@@ -159,6 +167,11 @@ public:
return implementation_->get_u((*eta_ast_prime_evaluations_)[entity_index]);
}
+ void get_u(const size_t entity_index, StateType& u) const
+ {
+ implementation_->get_u((*eta_ast_prime_evaluations_)[entity_index], u);
+ }
+
StateType get_alpha(const StateType& u) const
{
const auto alpha = implementation_->get_alpha(u)->first;
@@ -183,12 +196,25 @@ public:
densities_stencil, precomputed_fluxes, dd);
}
- void apply_inverse_hessian(const size_t entity_index, const StateType& u, StateType& Hinv_u) const
+ void apply_kinetic_flux_with_kinetic_reconstruction(const RangeFieldType& h_inv,
+ const QuadratureWeightsType& densities_left,
+ const QuadratureWeightsType& densities_entity,
+ const QuadratureWeightsType& densities_right,
+ StateType* u_left,
+ StateType* u_entity,
+ StateType* u_right,
+ const size_t dd) const
{
- implementation_->apply_inverse_hessian((*eta_ast_twoprime_evaluations_)[entity_index], u, Hinv_u);
+ implementation_->template apply_kinetic_flux_with_kinetic_reconstruction<slope>(
+ h_inv, densities_left, densities_entity, densities_right, u_left, u_entity, u_right, dd);
}
- void store_evaluations(const size_t entity_index, StateType& alpha, bool check = true)
+ void apply_inverse_hessian(const size_t entity_index, StateType& u) const
+ {
+ implementation_->apply_inverse_hessian((*eta_ast_twoprime_evaluations_)[entity_index], u);
+ }
+
+ void store_evaluations(const size_t entity_index, StateType& alpha, const RangeFieldType psi_min, bool check = true)
{
implementation_->store_exp_evaluations(exp_evaluations_[entity_index], alpha);
if (entropy != EntropyType::MaxwellBoltzmann) {
@@ -199,15 +225,15 @@ public:
set_eta_ast_pointers();
// check for inf and nan and very low densities
if (check) {
- const auto& u = get_u(entity_index);
- for (auto&& entry : u)
- if (std::isnan(entry) || std::isinf(entry))
- DUNE_THROW(Dune::MathError, "inf or nan in u!");
- const auto rho = basis_functions().density(u);
- if (rho < 1e-9) {
- alpha = basis_functions().alpha_iso(1e-8);
- store_evaluations(entity_index, alpha, false);
- }
+ const auto u = get_u(entity_index);
+ const double* u_ptr = &(u[0]);
+ // if (std::any_of(u_ptr, u_ptr + basis_dimRange, [](const auto& a) { return std::isnan(a) || std::isinf(a); }))
+ const auto val = XT::Common::reduce(u_ptr, u_ptr + basis_dimRange, 0.);
+ if (std::isnan(val) || std::isinf(val))
+ DUNE_THROW(Dune::MathError, "inf or nan in u!");
+ const bool changed = basis_functions().adjust_alpha_to_ensure_min_density(alpha, psi_min);
+ if (changed)
+ store_evaluations(entity_index, alpha, psi_min, false);
}
}
@@ -297,6 +323,11 @@ public:
return boundary_distribution_evaluations_;
}
+ RangeReturnType evaluate_with_alpha(const StateType& alpha) const
+ {
+ return implementation_->evaluate_with_alpha(alpha);
+ }
+
private:
std::shared_ptr<ImplementationType> implementation_;
std::vector<QuadratureWeightsType> exp_evaluations_;
diff --git a/dune/gdt/test/momentmodels/mn-discretization.hh b/dune/gdt/test/momentmodels/mn-discretization.hh
index 038cad0a2..6f8f9703e 100644
--- a/dune/gdt/test/momentmodels/mn-discretization.hh
+++ b/dune/gdt/test/momentmodels/mn-discretization.hh
@@ -148,14 +148,15 @@ struct HyperbolicMnDiscretization
using BoundaryOperator =
LocalAdvectionFvBoundaryTreatmentByCustomExtrapolationOperator<I, VectorType, GV, dimRange>;
using LambdaType = typename BoundaryOperator::LambdaType;
- using StateType = typename BoundaryOperator::StateType;
+ using DynamicStateType = typename BoundaryOperator::DynamicStateType;
LambdaType boundary_lambda =
[&boundary_values](const I& intersection,
const FieldVector<RangeFieldType, dimDomain - 1>& xx_in_reference_intersection_coordinates,
const AnalyticalFluxType& /*flux*/,
- const StateType& /*u*/,
- const XT::Common::Parameter& /*param*/) {
- return boundary_values->evaluate(intersection.geometry().global(xx_in_reference_intersection_coordinates));
+ const DynamicStateType& /*u*/,
+ DynamicStateType& v,
+ const XT::Common::Parameter& param) {
+ boundary_values->evaluate(intersection.geometry().global(xx_in_reference_intersection_coordinates), v, param);
};
XT::Grid::ApplyOn::NonPeriodicBoundaryIntersections<GV> filter;
advection_operator.append(boundary_lambda, {}, filter);
@@ -243,13 +244,13 @@ struct HyperbolicMnTest
void run(const double tol = TestCaseType::ExpectedResultsType::tol)
{
auto norms = HyperbolicMnDiscretization<TestCaseType>::run(
- 1,
+ DXTC_CONFIG.get("num_save_steps", 1),
0,
TestCaseType::quad_order,
TestCaseType::quad_refinements,
- "",
+ DXTC_CONFIG.get("grid_size", ""),
2,
- TestCaseType::t_end,
+ DXTC_CONFIG.get("t_end", TestCaseType::t_end),
"test",
Dune::GDT::is_full_moment_basis<typename TestCaseType::MomentBasis>::value)
.first;
diff --git a/dune/gdt/test/momentmodels/pn-discretization.hh b/dune/gdt/test/momentmodels/pn-discretization.hh
index cf2edc275..10c40bb69 100644
--- a/dune/gdt/test/momentmodels/pn-discretization.hh
+++ b/dune/gdt/test/momentmodels/pn-discretization.hh
@@ -206,8 +206,9 @@ struct HyperbolicPnDiscretization
using BoundaryValueType = typename ProblemType::BoundaryValueType;
static constexpr size_t dimDomain = MomentBasis::dimDomain;
static constexpr size_t dimRange = MomentBasis::dimRange;
- using MatrixType = typename XT::LA::Container<RangeFieldType>::MatrixType;
- using VectorType = typename XT::LA::Container<RangeFieldType>::VectorType;
+ static const auto la_backend = TestCaseType::la_backend;
+ using MatrixType = typename XT::LA::Container<RangeFieldType, la_backend>::MatrixType;
+ using VectorType = typename XT::LA::Container<RangeFieldType, la_backend>::VectorType;
//******************* create grid and FV space ***************************************
auto grid_config = ProblemType::default_grid_cfg();
@@ -294,14 +295,15 @@ struct HyperbolicPnDiscretization
using BoundaryOperator =
LocalAdvectionFvBoundaryTreatmentByCustomExtrapolationOperator<I, VectorType, GV, dimRange>;
using LambdaType = typename BoundaryOperator::LambdaType;
- using StateType = typename BoundaryOperator::StateType;
+ using DynamicStateType = typename BoundaryOperator::DynamicStateType;
LambdaType boundary_lambda =
[&boundary_values](const I& intersection,
const FieldVector<RangeFieldType, dimDomain - 1>& xx_in_reference_intersection_coordinates,
const AnalyticalFluxType& /*flux*/,
- const StateType& /*u*/,
- const XT::Common::Parameter& /*param*/) {
- return boundary_values->evaluate(intersection.geometry().global(xx_in_reference_intersection_coordinates));
+ const DynamicStateType& /*u*/,
+ DynamicStateType& v,
+ const XT::Common::Parameter& param) {
+ boundary_values->evaluate(intersection.geometry().global(xx_in_reference_intersection_coordinates), v, param);
};
XT::Grid::ApplyOn::NonPeriodicBoundaryIntersections<GV> filter;
advection_operator.append(boundary_lambda, {}, filter);
--
GitLab