diff --git a/dune/gdt/local/operators/advection-fv.hh b/dune/gdt/local/operators/advection-fv.hh
index fdebe9549ea15be3303b34725704aa7d06fb21fb..36f59a3a15b4056165ff7416d753b8034115fef2 100644
--- a/dune/gdt/local/operators/advection-fv.hh
+++ b/dune/gdt/local/operators/advection-fv.hh
@@ -201,10 +201,12 @@ public:
using typename BaseType::SourceSpaceType;
using typename BaseType::SourceType;
+ using D = typename IntersectionType::ctype;
using StateDomainType = FieldVector<typename SGV::ctype, SGV::dimension>;
using DynamicStateType = typename CouplingOperatorType::DynamicStateType;
- using LambdaType = std::function<void(const DynamicStateType& /*u*/,
- const StateDomainType& /*n*/,
+ using LambdaType = std::function<void(const IntersectionType& /*intersection*/,
+ const FieldVector<D, d - 1>& /*xx_in_reference_intersection_coordinates*/,
+ const DynamicStateType& /*u*/,
DynamicStateType& /*g*/,
const XT::Common::Parameter& /*param*/)>;
@@ -286,8 +288,8 @@ public:
local_sources_[0]->evaluate(source_is_elementwise_constant_ ? CouplingOperatorType::static_x
: intersection().geometryInInside().center(),
u_);
- const auto normal = intersection().centerUnitOuterNormal();
- numerical_boundary_flux_(u_, normal, g_, param);
+ numerical_boundary_flux_(
+ intersection(), intersection().geometry().local(intersection().geometry().center()), u_, 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)
diff --git a/dune/gdt/test/momentmodels/entropic-coords-mn-discretization.hh b/dune/gdt/test/momentmodels/entropic-coords-mn-discretization.hh
index 4b03770ebaea7ad52e4b10762f46e5f520022b4f..9a76b1beab85844055f0d1411c077ba081ae3fe8 100644
--- a/dune/gdt/test/momentmodels/entropic-coords-mn-discretization.hh
+++ b/dune/gdt/test/momentmodels/entropic-coords-mn-discretization.hh
@@ -137,13 +137,15 @@ struct HyperbolicEntropicCoordsMnDiscretization
const auto u_local_func = u.local_discrete_function();
const auto alpha_local_func = alpha.local_discrete_function();
+ const auto entropy_local_func = entropy_flux->derived_local_function();
XT::Common::FieldVector<RangeFieldType, dimRange> u_local;
for (auto&& element : Dune::elements(grid_view)) {
u_local_func->bind(element);
alpha_local_func->bind(element);
+ entropy_local_func->bind(element);
for (size_t ii = 0; ii < dimRange; ++ii)
u_local[ii] = u_local_func->dofs().get_entry(ii);
- const auto alpha_local = analytical_flux->get_alpha(u_local);
+ const auto alpha_local = entropy_local_func->get_alpha(u_local, false)->first;
for (size_t ii = 0; ii < dimRange; ++ii)
alpha_local_func->dofs().set_entry(ii, alpha_local[ii]);
}
@@ -309,7 +311,7 @@ struct HyperbolicEntropicCoordsMnDiscretization
// The hessian has entries in the order of psi_min, the inverse thus scales with 1/psi_min, and thus the timestep
// should be psi_min to get an update of order 1
- double initial_dt = 1e-14; // std::min(dt, min_acceptable_density);
+ double initial_dt = dx / 100.; // std::min(dt, min_acceptable_density);
timestepper.solve(t_end,
initial_dt,
num_save_steps,
diff --git a/dune/gdt/test/momentmodels/entropyflux.hh b/dune/gdt/test/momentmodels/entropyflux.hh
index bbd3e6f7e8372ffff3a5aeb5e761b7b233ee54d0..df43f947c7047a9798706f924f7f3d1c1f829529 100644
--- a/dune/gdt/test/momentmodels/entropyflux.hh
+++ b/dune/gdt/test/momentmodels/entropyflux.hh
@@ -331,7 +331,6 @@ public:
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
@@ -342,6 +341,21 @@ public:
return implementation_->get_alpha(u, *alpha_iso, regularize);
}
+ std::unique_ptr<AlphaReturnType> get_alpha(const E& entity, const StateType& u, const bool regularize) const
+ {
+ thread_local const auto local_func = derived_local_function();
+ local_func->bind(entity);
+ return local_func->get_alpha(u, regularize);
+ }
+
+ StateType evaluate_kinetic_outflow(const typename AlphaReturnType::first_type& alpha_i,
+ const DomainType& n_ij,
+ const size_t dd) const
+
+ {
+ return implementation_->evaluate_kinetic_outflow(alpha_i, n_ij, dd);
+ } // DomainType evaluate_kinetic_outflow(...)
+
const MomentBasis& basis_functions() const
{
return implementation_->basis_functions();
diff --git a/dune/gdt/test/momentmodels/entropyflux_implementations.hh b/dune/gdt/test/momentmodels/entropyflux_implementations.hh
index fc55f30a87a3773590d957cf9ca038822ea06df0..060e428eca2fc396f5feb68eac435242f7a90081 100644
--- a/dune/gdt/test/momentmodels/entropyflux_implementations.hh
+++ b/dune/gdt/test/momentmodels/entropyflux_implementations.hh
@@ -545,6 +545,25 @@ public:
return ret;
} // DomainType evaluate_kinetic_flux_with_alphas(...)
+ DomainType evaluate_kinetic_outflow(const DomainType& alpha_i, const FluxDomainType& n_ij, const size_t dd) const
+
+ {
+ auto& eta_ast_prime_vals_i = working_storage();
+ evaluate_eta_ast_prime(alpha_i, M_, eta_ast_prime_vals_i);
+ DomainType ret(0);
+ for (size_t ll = 0; ll < quad_points_.size(); ++ll) {
+ const auto position = quad_points_[ll][dd];
+ if (position * n_ij[dd] > 0.) {
+ RangeFieldType factor = eta_ast_prime_vals_i[ll] * quad_weights_[ll] * position;
+ const auto* basis_ll = &(M_.get_entry_ref(ll, 0.));
+ for (size_t ii = 0; ii < basis_dimRange; ++ii)
+ ret[ii] += basis_ll[ii] * factor;
+ }
+ } // ll
+ ret *= n_ij[dd];
+ return ret;
+ } // DomainType evaluate_kinetic_outflow(...)
+
// Calculates left and right kinetic flux with reconstructed densities. Ansatz distribution values contains
// evaluations of the ansatz distribution at each quadrature point for a stencil of three entities. The distributions
// are reconstructed pointwise for each quadrature point and the resulting (part of) the kinetic flux is <
@@ -2051,6 +2070,46 @@ public:
return ret;
} // DomainType evaluate_kinetic_flux(...)
+ // DomainType evaluate_kinetic_outflow(const DomainType& alpha,
+ // const FluxDomainType& n_ij,
+ // const size_t dd) const
+ // {
+ // return evaluate_kinetic_outflow(alpha, n_ij, dd);
+ // } // DomainType evaluate_kinetic_outflow(...)
+
+ DomainType evaluate_kinetic_outflow(const BlockVectorType& alpha_i, const FluxDomainType& n_ij, const size_t dd) const
+ {
+ // calculate \sum_{i=1}^d < \omega_i m G_\alpha(u) > n_i
+ auto& eta_ast_prime_vals = working_storage();
+ evaluate_eta_ast_prime(alpha_i, M_, eta_ast_prime_vals);
+ DomainType ret(0);
+ for (size_t jj = 0; jj < num_blocks; ++jj) {
+ const auto offset = block_size * jj;
+ if (quad_signs_[jj][dd] == 0) {
+ // in this case, we have to decide which eta_ast_prime_vals to use for each quadrature point individually
+ for (size_t ll = 0; ll < quad_points_[jj].size(); ++ll) {
+ const auto position = quad_points_[jj][ll][dd];
+ if (position * n_ij[dd] > 0.) {
+ RangeFieldType factor = eta_ast_prime_vals[jj][ll] * quad_weights_[jj][ll] * position;
+ for (size_t ii = 0; ii < block_size; ++ii)
+ ret[offset + ii] += M_[jj].get_entry(ll, ii) * factor;
+ }
+ } // ll
+ } else {
+ // all quadrature points have the same sign
+ if (quad_signs_[jj][dd] * n_ij[dd] > 0.) {
+ for (size_t ll = 0; ll < quad_points_[jj].size(); ++ll) {
+ const RangeFieldType factor = eta_ast_prime_vals[jj][ll] * quad_weights_[jj][ll] * quad_points_[jj][ll][dd];
+ for (size_t ii = 0; ii < block_size; ++ii)
+ ret[offset + ii] += M_[jj].get_entry(ll, ii) * factor;
+ } // ll
+ }
+ } // quad_sign
+ } // jj
+ ret *= n_ij[dd];
+ return ret;
+ } // DomainType evaluate_kinetic_outflow(...)
+
// Calculates left and right kinetic flux with reconstructed densities. Ansatz distribution values contains
// evaluations of the ansatz distribution at each quadrature point for a stencil of three entities. The distributions
// are reconstructed pointwise for each quadrature point and the resulting (part of) the kinetic flux is <
@@ -3086,6 +3145,37 @@ public:
return ret;
} // DomainType evaluate_kinetic_flux(...)
+ DomainType evaluate_kinetic_outflow(const DomainType& alpha_i, const FluxDomainType& n_ij, const size_t dd) const
+ {
+ return evaluate_kinetic_outflow(XT::LA::convert_to<VectorType>(alpha_i), n_ij, dd);
+ }
+
+ DomainType evaluate_kinetic_outflow(const VectorType& alpha_i, const FluxDomainType& n_ij, const size_t dd) const
+ {
+ auto& eta_ast_prime_vals = working_storage();
+ 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 auto& faces = basis_functions_.triangulation().faces();
+ LocalVectorType local_ret;
+ for (size_t jj = 0; jj < num_faces_; ++jj) {
+ const bool positive_dir = v_positive_[jj][dd];
+ if ((n_ij[dd] > 0. && positive_dir) || (n_ij[dd] < 0. && !positive_dir)) {
+ local_ret *= 0.;
+ const auto& vertices = faces[jj]->vertices();
+ for (size_t ll = 0; ll < quad_weights_[jj].size(); ++ll) {
+ RangeFieldType factor = eta_ast_prime_vals[jj][ll] * quad_weights_[jj][ll] * quad_points_[jj][ll][dd];
+ for (size_t ii = 0; ii < 3; ++ii)
+ local_ret[ii] += M_[jj][ll][ii] * factor;
+ } // ll (quad points)
+ for (size_t ii = 0; ii < 3; ++ii)
+ ret[vertices[ii]->index()] += local_ret[ii];
+ }
+ } // jj (faces)
+ ret *= n_ij[dd];
+ return ret;
+ } // DomainType evaluate_kinetic_outflow(...)
+
// Calculates left and right kinetic flux with reconstructed densities. Ansatz distribution values contains
// evaluations of the ansatz distribution at each quadrature point for a stencil of three entities. The distributions
// are reconstructed pointwise for each quadrature point and the resulting (part of) the kinetic flux is <
@@ -4474,6 +4564,26 @@ public:
return ret;
} // DomainType evaluate_kinetic_flux(...)
+ DomainType evaluate_kinetic_outflow(const VectorType& alpha_i, const FluxDomainType& n_ij, const size_t dd) const
+ {
+ auto& eta_ast_prime_vals = working_storage();
+ 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;
+ if (n_ij[dd] < 0.) {
+ for (size_t ll = 0; ll < first_positive; ++ll)
+ ret[ll] = eta_ast_prime_vals[ll] * grid_points_[ll];
+ } else {
+ for (size_t ll = first_positive; ll < dimRange; ++ll)
+ ret[ll] = second_exp_vals[ll] * grid_points_[ll];
+ }
+ ret *= n_ij[dd] * interval_length;
+ ret[0] *= 0.5;
+ ret[dimRange - 1] *= 0.5;
+ return ret;
+ } // DomainType evaluate_kinetic_outflow(...)
+
// Calculates left and right kinetic flux with reconstructed densities. Ansatz distribution values contains
// evaluations of the ansatz distribution at each quadrature point for a stencil of three entities. The distributions
// are reconstructed pointwise for each quadrature point and the resulting (part of) the kinetic flux is
@@ -5178,6 +5288,27 @@ public:
return ret;
} // DomainType evaluate_kinetic_flux(...)
+ DomainType evaluate_kinetic_outflow(const VectorType& alpha_i, const FluxDomainType& n_ij, const size_t dd) const
+ {
+ assert(dd == 0);
+ auto& eta_ast_prime_vals = working_storage();
+ 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);
+ for (size_t jj = 0; jj < num_intervals; ++jj) {
+ for (size_t ll = 0; ll < quad_weights_[jj].size(); ++ll) {
+ const auto position = quad_points_[jj][ll];
+ if (position * n_ij[dd] > 0.) {
+ const RangeFieldType factor = eta_ast_prime_vals[jj][ll] * quad_weights_[jj][ll] * position;
+ for (size_t ii = 0; ii < 2; ++ii)
+ ret[jj + ii] += M_[jj][ll][ii] * factor;
+ }
+ } // ll (quad points)
+ } // jj (faces)
+ ret *= n_ij[dd];
+ return ret;
+ } // DomainType evaluate_kinetic_flux(...)
+
// Calculates left and right kinetic flux with reconstructed densities. Ansatz distribution values contains
// evaluations of the ansatz distribution at each quadrature point for a stencil of three entities. The distributions
// are reconstructed pointwise for each quadrature point and the resulting (part of) the kinetic flux is <
diff --git a/dune/gdt/test/momentmodels/mn-discretization.hh b/dune/gdt/test/momentmodels/mn-discretization.hh
index 9f78fde32142e040f8ade3f0a92e20e1fd6adc30..634f1691b41ebb42e2a3f51043339f9af09d73c6 100644
--- a/dune/gdt/test/momentmodels/mn-discretization.hh
+++ b/dune/gdt/test/momentmodels/mn-discretization.hh
@@ -66,8 +66,10 @@ struct HyperbolicMnDiscretization
static constexpr size_t dimDomain = MomentBasis::dimDomain;
static constexpr size_t dimRange = MomentBasis::dimRange;
static const auto la_backend = TestCaseType::la_backend;
+ using DomainType = FieldVector<RangeFieldType, dimDomain>;
using MatrixType = typename XT::LA::Container<RangeFieldType, la_backend>::MatrixType;
using VectorType = typename XT::LA::Container<RangeFieldType, la_backend>::VectorType;
+ using DynamicRangeType = DynamicVector<RangeFieldType>;
//******************* create grid and FV space ***************************************
auto grid_config = ProblemType::default_grid_cfg();
@@ -94,6 +96,7 @@ struct HyperbolicMnDiscretization
using AnalyticalFluxType = typename ProblemType::FluxType;
using EntropyFluxType = typename ProblemType::ActualFluxType;
auto analytical_flux = problem.flux();
+ auto* entropy_flux = dynamic_cast<EntropyFluxType*>(analytical_flux.get());
// 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.
@@ -146,24 +149,60 @@ struct HyperbolicMnDiscretization
// boundary treatment
using BoundaryOperator =
- LocalAdvectionFvBoundaryTreatmentByCustomExtrapolationOperator<I, VectorType, GV, dimRange>;
+ // LocalAdvectionFvBoundaryTreatmentByCustomExtrapolationOperator<I, VectorType, GV, dimRange>;
+ LocalAdvectionFvBoundaryTreatmentByCustomNumericalFluxOperator<I, VectorType, GV, dimRange>;
using LambdaType = typename BoundaryOperator::LambdaType;
- using DynamicStateType = typename BoundaryOperator::DynamicStateType;
+ // 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
+ // DynamicStateType& /*u*/, DynamicStateType& v, const XT::Common::Parameter& param) {
+ // boundary_values->evaluate(intersection.geometry().global(xx_in_reference_intersection_coordinates), v,
+ // param);
+ // };
+
+ // store boundary influx 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 DynamicRangeType boundary_flux =
+ problem.kinetic_boundary_flux(x, intersection.centerUnitOuterNormal()[dd], dd);
+ boundary_fluxes.insert(std::make_pair(x, boundary_flux));
+ }
+ using GenericFunctionType = XT::Functions::GenericFunction<dimDomain, dimRange, 1, RangeFieldType>;
+ GenericFunctionType boundary_kinetic_fluxes(
+ 1, [&](const DomainType& x, DynamicRangeType& ret, const XT::Common::Parameter&) { ret = boundary_fluxes[x]; });
+
LambdaType boundary_lambda =
- [&boundary_values](const I& intersection,
- const FieldVector<RangeFieldType, dimDomain - 1>& xx_in_reference_intersection_coordinates,
- const AnalyticalFluxType& /*flux*/,
- const DynamicStateType& /*u*/,
- DynamicStateType& v,
- const XT::Common::Parameter& param) {
- boundary_values->evaluate(intersection.geometry().global(xx_in_reference_intersection_coordinates), v, param);
+ [&boundary_kinetic_fluxes,
+ &entropy_flux](const I& intersection,
+ const FieldVector<RangeFieldType, dimDomain - 1>& xx_in_reference_intersection_coordinates,
+ const DynamicRangeType& u,
+ DynamicRangeType& g,
+ const XT::Common::Parameter& param) {
+ // influx
+ boundary_kinetic_fluxes.evaluate(
+ intersection.geometry().global(xx_in_reference_intersection_coordinates), g, param);
+ // outflux
+ const auto& entity = intersection.inside();
+ const auto dd = intersection.indexInInside() / 2;
+ const auto alpha_entity = entropy_flux->get_alpha(entity, u, true)->first;
+ const auto outflux =
+ entropy_flux->evaluate_kinetic_outflow(alpha_entity, intersection.centerUnitOuterNormal(), dd);
+ g -= outflux;
+ g *= -1;
};
+
XT::Grid::ApplyOn::NonPeriodicBoundaryIntersections<GV> filter;
advection_operator.append(boundary_lambda, {}, filter);
constexpr double epsilon = 1e-11;
auto slope = TestCaseType::RealizabilityLimiterChooserType::template make_slope<EigenvectorWrapperType>(
- *dynamic_cast<EntropyFluxType*>(analytical_flux.get()), *basis_functions, epsilon);
+ *entropy_flux, *basis_functions, epsilon);
ReconstructionOperatorType reconstruction_operator(*analytical_flux, *boundary_values, fv_space, *slope, false);
ReconstructionAdvectionOperatorType reconstruction_advection_operator(advection_operator, reconstruction_operator);