diff --git a/dune/gdt/local/operators/advection-dg.hh b/dune/gdt/local/operators/advection-dg.hh
index 49b42129ac585e1087b52797491a0cea4e3efe4a..a978ecbaf5be2111280b9b366eb7af35e9fbab73 100644
--- a/dune/gdt/local/operators/advection-dg.hh
+++ b/dune/gdt/local/operators/advection-dg.hh
@@ -15,6 +15,7 @@
#include <dune/geometry/quadraturerules.hh>
+#include <dune/xt/common/memory.hh>
#include <dune/xt/common/parameter.hh>
#include <dune/gdt/exceptions.hh>
@@ -22,6 +23,7 @@
#include <dune/gdt/type_traits.hh>
#include <dune/gdt/local/numerical-fluxes/interface.hh>
+#include <dune/gdt/tools/local-mass-matrix.hh>
#include "interfaces.hh"
@@ -47,15 +49,24 @@ public:
using typename BaseType::SourceType;
using FluxType = XT::Functions::FunctionInterface<m, d, m, RF>;
+ using LocalMassMatrixProviderType = LocalMassMatrixProvider<RGV, m, 1, RF>;
LocalAdvectionDgVolumeOperator(const FluxType& flux)
: BaseType(flux.parameter_type())
, flux_(flux)
{}
+ /// Applies the inverse of the local mass matrix.
+ LocalAdvectionDgVolumeOperator(const LocalMassMatrixProviderType& local_mass_matrices, const FluxType& flux)
+ : BaseType(flux.parameter_type())
+ , flux_(flux)
+ , local_mass_matrices_(local_mass_matrices)
+ {}
+
LocalAdvectionDgVolumeOperator(const ThisType& other)
: BaseType(other.parameter_type())
, flux_(other.flux_)
+ , local_mass_matrices_(other.local_mass_matrices_)
{}
std::unique_ptr<BaseType> copy() const override final
@@ -69,6 +80,8 @@ public:
{
const auto& element = local_range.element();
const auto& basis = local_range.basis();
+ local_dofs_.resize(basis.size(param));
+ local_dofs_ *= 0.;
const auto local_source = source.local_discrete_function(element);
const auto local_source_order = local_source->order(param);
const auto local_basis_order = basis.order(param);
@@ -84,13 +97,21 @@ public:
const auto flux_value = flux_.evaluate(source_value, param);
// compute
for (size_t ii = 0; ii < basis.size(param); ++ii)
- local_range.dofs()[ii] += integration_factor * quadrature_weight * -1. * (flux_value * basis_jacobians_[ii]);
+ local_dofs_[ii] += integration_factor * quadrature_weight * -1. * (flux_value * basis_jacobians_[ii]);
}
+ // apply local mass matrix, if required (not optimal, uses a temporary)
+ if (local_mass_matrices_.valid())
+ local_dofs_ = local_mass_matrices_.access().local_mass_matrix_inverse(element) * local_dofs_;
+ // add to local range
+ for (size_t ii = 0; ii < basis.size(param); ++ii)
+ local_range.dofs()[ii] += local_dofs_[ii];
} // ... apply(...)
private:
const FluxType& flux_;
+ const XT::Common::ConstStorageProvider<LocalMassMatrixProviderType> local_mass_matrices_;
mutable std::vector<typename LocalRangeType::LocalBasisType::DerivativeRangeType> basis_jacobians_;
+ mutable XT::LA::CommonDenseVector<RF> local_dofs_;
}; // class LocalAdvectionDgVolumeOperator
@@ -104,17 +125,17 @@ template <class I,
class SGV,
size_t m = 1,
class SR = double,
- class RR = SR,
+ class IRR = SR,
class IRGV = SGV,
class IRV = SV,
- class ORR = RR,
+ class ORR = IRR,
class ORGV = IRGV,
class ORV = IRV>
class LocalAdvectionDgCouplingOperator
- : public LocalIntersectionOperatorInterface<I, SV, SGV, m, 1, SR, m, 1, RR, IRGV, IRV, ORGV, ORV>
+ : public LocalIntersectionOperatorInterface<I, SV, SGV, m, 1, SR, m, 1, IRR, IRGV, IRV, ORGV, ORV>
{
- using ThisType = LocalAdvectionDgCouplingOperator<I, SV, SGV, m, SR, RR, IRGV, IRV, ORR, ORGV, ORV>;
- using BaseType = LocalIntersectionOperatorInterface<I, SV, SGV, m, 1, SR, m, 1, RR, IRGV, IRV, ORGV, ORV>;
+ using ThisType = LocalAdvectionDgCouplingOperator<I, SV, SGV, m, SR, IRR, IRGV, IRV, ORR, ORGV, ORV>;
+ using BaseType = LocalIntersectionOperatorInterface<I, SV, SGV, m, 1, SR, m, 1, IRR, IRGV, IRV, ORGV, ORV>;
public:
using BaseType::d;
@@ -124,7 +145,8 @@ public:
using typename BaseType::LocalOutsideRangeType;
using typename BaseType::SourceType;
- using NumericalFluxType = NumericalFluxInterface<d, m, RR>;
+ using NumericalFluxType = NumericalFluxInterface<d, m, IRR>;
+ using LocalMassMatrixProviderType = LocalMassMatrixProvider<IRGV, m, 1, IRR>;
LocalAdvectionDgCouplingOperator(const NumericalFluxType& numerical_flux, bool compute_outside = true)
: BaseType(numerical_flux.parameter_type())
@@ -132,10 +154,21 @@ public:
, compute_outside_(compute_outside)
{}
+ /// Applies the inverse of the local mass matrix.
+ LocalAdvectionDgCouplingOperator(const LocalMassMatrixProviderType& local_mass_matrices,
+ const NumericalFluxType& numerical_flux,
+ bool compute_outside = true)
+ : BaseType(numerical_flux.parameter_type())
+ , numerical_flux_(numerical_flux.copy())
+ , compute_outside_(compute_outside)
+ , local_mass_matrices_(local_mass_matrices)
+ {}
+
LocalAdvectionDgCouplingOperator(const ThisType& other)
: BaseType(other.parameter_type())
, numerical_flux_(other.numerical_flux_->copy())
, compute_outside_(other.compute_outside_)
+ , local_mass_matrices_(other.local_mass_matrices_)
{}
std::unique_ptr<BaseType> copy() const override final
@@ -153,6 +186,10 @@ public:
const auto& outside_element = local_range_outside.element();
const auto& inside_basis = local_range_inside.basis();
const auto& outside_basis = local_range_outside.basis();
+ inside_local_dofs_.resize(inside_basis.size(param));
+ outside_local_dofs_.resize(outside_basis.size(param));
+ inside_local_dofs_ *= 0.;
+ outside_local_dofs_ *= 0.;
const auto u = source.local_discrete_function(inside_element);
const auto v = source.local_discrete_function(outside_element);
const auto integrand_order = std::max(inside_basis.order(param), outside_basis.order(param))
@@ -178,18 +215,33 @@ public:
param);
// compute
for (size_t ii = 0; ii < inside_basis.size(param); ++ii)
- local_range_inside.dofs()[ii] += integration_factor * quadrature_weight * (g * inside_basis_values_[ii]);
+ inside_local_dofs_[ii] += integration_factor * quadrature_weight * (g * inside_basis_values_[ii]);
if (compute_outside_)
for (size_t ii = 0; ii < outside_basis.size(param); ++ii)
- local_range_outside.dofs()[ii] -= integration_factor * quadrature_weight * (g * outside_basis_values_[ii]);
+ outside_local_dofs_[ii] -= integration_factor * quadrature_weight * (g * outside_basis_values_[ii]);
}
+ // apply local mass matrix, if required (not optimal, uses a temporary)
+ if (local_mass_matrices_.valid())
+ inside_local_dofs_ = local_mass_matrices_.access().local_mass_matrix_inverse(inside_element) * inside_local_dofs_;
+ if (compute_outside_ && local_mass_matrices_.valid())
+ outside_local_dofs_ =
+ local_mass_matrices_.access().local_mass_matrix_inverse(outside_element) * outside_local_dofs_;
+ // add to local range
+ for (size_t ii = 0; ii < inside_basis.size(param); ++ii)
+ local_range_inside.dofs()[ii] += inside_local_dofs_[ii];
+ if (compute_outside_)
+ for (size_t ii = 0; ii < outside_basis.size(param); ++ii)
+ local_range_outside.dofs()[ii] += outside_local_dofs_[ii];
} // ... apply(...)
private:
const std::unique_ptr<const NumericalFluxType> numerical_flux_;
const bool compute_outside_;
+ const XT::Common::ConstStorageProvider<LocalMassMatrixProviderType> local_mass_matrices_;
mutable std::vector<typename LocalInsideRangeType::LocalBasisType::RangeType> inside_basis_values_;
mutable std::vector<typename LocalOutsideRangeType::LocalBasisType::RangeType> outside_basis_values_;
+ mutable XT::LA::CommonDenseVector<IRR> inside_local_dofs_;
+ mutable XT::LA::CommonDenseVector<ORR> outside_local_dofs_;
}; // class LocalAdvectionDgCouplingOperator
@@ -217,6 +269,7 @@ public:
using StateRangeType = typename XT::Functions::RangeTypeSelector<SF, m, 1>::type;
using LambdaType = std::function<StateRangeType(
const StateRangeType& /*u*/, const StateDomainType& /*n*/, const XT::Common::Parameter& /*param*/)>;
+ using LocalMassMatrixProviderType = LocalMassMatrixProvider<RGV, m, 1, RF>;
LocalAdvectionDgBoundaryTreatmentByCustomNumericalFluxOperator(
LambdaType numerical_boundary_flux,
@@ -227,10 +280,23 @@ public:
, numerical_flux_order_(numerical_flux_order)
{}
+ /// Applies the inverse of the local mass matrix.
+ LocalAdvectionDgBoundaryTreatmentByCustomNumericalFluxOperator(
+ const LocalMassMatrixProviderType& local_mass_matrices,
+ LambdaType numerical_boundary_flux,
+ const int numerical_flux_order,
+ const XT::Common::ParameterType& numerical_flux_param_type = {})
+ : BaseType(numerical_flux_param_type)
+ , numerical_boundary_flux_(numerical_boundary_flux)
+ , numerical_flux_order_(numerical_flux_order)
+ , local_mass_matrices_(local_mass_matrices)
+ {}
+
LocalAdvectionDgBoundaryTreatmentByCustomNumericalFluxOperator(const ThisType& other)
: BaseType(other.parameter_type())
, numerical_boundary_flux_(other.numerical_boundary_flux_)
, numerical_flux_order_(other.numerical_flux_order_)
+ , local_mass_matrices_(other.local_mass_matrices_)
{}
std::unique_ptr<BaseType> copy() const override final
@@ -246,6 +312,8 @@ public:
{
const auto& element = local_range_inside.element();
const auto& inside_basis = local_range_inside.basis();
+ inside_local_dofs_.resize(inside_basis.size(param));
+ inside_local_dofs_ *= 0.;
const auto u = source.local_discrete_function(element);
const auto integrand_order = inside_basis.order(param) + numerical_flux_order_ * u->order(param);
for (const auto& quadrature_point :
@@ -262,14 +330,22 @@ public:
const auto g = numerical_boundary_flux_(u->evaluate(point_in_inside_reference_element), normal, param);
// compute
for (size_t ii = 0; ii < inside_basis.size(param); ++ii)
- local_range_inside.dofs()[ii] += integration_factor * quadrature_weight * (g * inside_basis_values_[ii]);
+ inside_local_dofs_[ii] += integration_factor * quadrature_weight * (g * inside_basis_values_[ii]);
}
+ // apply local mass matrix, if required (not optimal, uses a temporary)
+ if (local_mass_matrices_.valid())
+ inside_local_dofs_ = local_mass_matrices_.access().local_mass_matrix_inverse(element) * inside_local_dofs_;
+ // add to local range
+ for (size_t ii = 0; ii < inside_basis.size(param); ++ii)
+ local_range_inside.dofs()[ii] += inside_local_dofs_[ii];
} // ... apply(...)
private:
const LambdaType numerical_boundary_flux_;
const int numerical_flux_order_;
+ const XT::Common::ConstStorageProvider<LocalMassMatrixProviderType> local_mass_matrices_;
mutable std::vector<typename LocalInsideRangeType::LocalBasisType::RangeType> inside_basis_values_;
+ mutable XT::LA::CommonDenseVector<RF> inside_local_dofs_;
}; // class LocalAdvectionDgBoundaryTreatmentByCustomNumericalFluxOperator
@@ -302,20 +378,34 @@ public:
const FluxType& /*flux*/,
const StateRangeType& /*u*/,
const XT::Common::Parameter& /*param*/)>;
+ using LocalMassMatrixProviderType = LocalMassMatrixProvider<RGV, m, 1, RF>;
+
+ LocalAdvectionDgBoundaryTreatmentByCustomExtrapolationOperator(
+ const NumericalFluxType& numerical_flux,
+ LambdaType boundary_extrapolation_lambda,
+ const XT::Common::ParameterType& boundary_treatment_param_type = {})
+ : BaseType(numerical_flux.parameter_type() + boundary_treatment_param_type)
+ , numerical_flux_(numerical_flux.copy())
+ , extrapolate_(boundary_extrapolation_lambda)
+ {}
+ /// Applies the inverse of the local mass matrix.
LocalAdvectionDgBoundaryTreatmentByCustomExtrapolationOperator(
+ const LocalMassMatrixProviderType& local_mass_matrices,
const NumericalFluxType& numerical_flux,
LambdaType boundary_extrapolation_lambda,
const XT::Common::ParameterType& boundary_treatment_param_type = {})
: BaseType(numerical_flux.parameter_type() + boundary_treatment_param_type)
, numerical_flux_(numerical_flux.copy())
, extrapolate_(boundary_extrapolation_lambda)
+ , local_mass_matrices_(local_mass_matrices)
{}
LocalAdvectionDgBoundaryTreatmentByCustomExtrapolationOperator(const ThisType& other)
: BaseType(other.parameter_type())
, numerical_flux_(other.numerical_flux_->copy())
, extrapolate_(other.extrapolate_)
+ , local_mass_matrices_(other.local_mass_matrices_)
{}
std::unique_ptr<BaseType> copy() const override final
@@ -331,6 +421,8 @@ public:
{
const auto& element = local_range_inside.element();
const auto& inside_basis = local_range_inside.basis();
+ inside_local_dofs_.resize(inside_basis.size(param));
+ inside_local_dofs_ *= 0.;
const auto local_source = source.local_discrete_function(element);
const auto integrand_order =
inside_basis.order(param) + numerical_flux_->flux().order(param) * local_source->order(param);
@@ -350,14 +442,22 @@ public:
const auto g = numerical_flux_->apply(u, v, normal, param);
// compute
for (size_t ii = 0; ii < inside_basis.size(param); ++ii)
- local_range_inside.dofs()[ii] += integration_factor * quadrature_weight * (g * inside_basis_values_[ii]);
+ inside_local_dofs_[ii] += integration_factor * quadrature_weight * (g * inside_basis_values_[ii]);
}
+ // apply local mass matrix, if required (not optimal, uses a temporary)
+ if (local_mass_matrices_.valid())
+ inside_local_dofs_ = local_mass_matrices_.access().local_mass_matrix_inverse(element) * inside_local_dofs_;
+ // add to local range
+ for (size_t ii = 0; ii < inside_basis.size(param); ++ii)
+ local_range_inside.dofs()[ii] += inside_local_dofs_[ii];
} // ... apply(...)
private:
const std::unique_ptr<const NumericalFluxType> numerical_flux_;
const LambdaType extrapolate_;
+ const XT::Common::ConstStorageProvider<LocalMassMatrixProviderType> local_mass_matrices_;
mutable std::vector<typename LocalInsideRangeType::LocalBasisType::RangeType> inside_basis_values_;
+ mutable XT::LA::CommonDenseVector<RF> inside_local_dofs_;
}; // class LocalAdvectionDgBoundaryTreatmentByCustomExtrapolationOperator
diff --git a/dune/gdt/operators/advection-dg.hh b/dune/gdt/operators/advection-dg.hh
index 96a5193b09e7c57f4d35dfd1abd9f0894e76a618..7947ca3d2baadce054cb70e0d6f23d93a65b0c54 100644
--- a/dune/gdt/operators/advection-dg.hh
+++ b/dune/gdt/operators/advection-dg.hh
@@ -21,6 +21,7 @@
#include <dune/xt/grid/intersection.hh>
#include <dune/xt/grid/type_traits.hh>
#include <dune/xt/grid/filters.hh>
+#include <dune/xt/grid/walker.hh>
#include <dune/gdt/local/bilinear-forms/integrals.hh>
#include <dune/gdt/local/functionals/integrals.hh>
@@ -28,6 +29,7 @@
#include <dune/gdt/local/integrands/product.hh>
#include <dune/gdt/local/operators/advection-dg.hh>
#include <dune/gdt/spaces/l2/finite-volume.hh>
+#include <dune/gdt/tools/local-mass-matrix.hh>
#include "interfaces.hh"
#include "localizable-operator.hh"
@@ -80,7 +82,13 @@ public:
, source_space_(source_space)
, range_space_(range_space)
, periodicity_exception_(periodicity_exception.copy())
- {}
+ , local_mass_matrix_provider_(assembly_grid_view_, range_space_)
+ {
+ // we assemble these once, to be used in each apply later on
+ auto walker = XT::Grid::make_walker(assembly_grid_view_);
+ walker.append(local_mass_matrix_provider_);
+ walker.walk(/*use_tbb=*/true);
+ }
AdvectionDgOperator(ThisType&& source) = default;
@@ -109,7 +117,8 @@ public:
const XT::Grid::IntersectionFilter<SGV>& filter = XT::Grid::ApplyOn::BoundaryIntersections<SGV>())
{
boundary_treatments_by_custom_numerical_flux_.emplace_back(
- new BoundaryTreatmentByCustomNumericalFluxOperatorType(numerical_boundary_treatment_flux,
+ new BoundaryTreatmentByCustomNumericalFluxOperatorType(local_mass_matrix_provider_,
+ numerical_boundary_treatment_flux,
numerical_boundary_treatment_flux_order,
boundary_treatment_parameter_type),
filter.copy());
@@ -122,27 +131,39 @@ public:
{
boundary_treatments_by_custom_extrapolation_.emplace_back(
new BoundaryTreatmentByCustomExtrapolationOperatorType(
- *numerical_flux_, extrapolation, extrapolation_parameter_type),
+ local_mass_matrix_provider_, *numerical_flux_, extrapolation, extrapolation_parameter_type),
filter.copy());
return *this;
}
/// \}
-protected:
- void append_standard_contributions(LocalizableOperatorBase<SGV, V, m, 1, F, SGV, m, 1, F, RGV, V>& localizable_op,
- const XT::Common::Parameter& param) const
+ void apply(const VectorType& source, VectorType& range, const XT::Common::Parameter& param = {}) const override
{
+ // some checks
+ DUNE_THROW_IF(!source.valid(), Exceptions::operator_error, "source contains inf or nan!");
+ DUNE_THROW_IF(!(this->parameter_type() <= param.type()),
+ Exceptions::operator_error,
+ "this->parameter_type() = " << this->parameter_type() << "\n param.type() = " << param.type());
+ range.set_all(0);
+ auto source_function = make_discrete_function(this->source_space_, source);
+ auto range_function = make_discrete_function(this->range_space_, range);
+ // set up the actual operator
+ auto localizable_op = make_localizable_operator(this->assembly_grid_view_, source_function, range_function);
// element contributions
- localizable_op.append(LocalAdvectionDgVolumeOperator<V, SGV, m, F, F, RGV, V>(numerical_flux_->flux()), param);
+ localizable_op.append(
+ LocalAdvectionDgVolumeOperator<V, SGV, m, F, F, RGV, V>(local_mass_matrix_provider_, numerical_flux_->flux()),
+ param);
// contributions from inner intersections
- localizable_op.append(LocalAdvectionDgCouplingOperator<I, V, SGV, m, F, F, RGV, V>(*numerical_flux_),
+ localizable_op.append(LocalAdvectionDgCouplingOperator<I, V, SGV, m, F, F, RGV, V>(
+ local_mass_matrix_provider_, *numerical_flux_, /*compute_outside=*/false),
param,
- XT::Grid::ApplyOn::InnerIntersectionsOnce<SGV>());
+ XT::Grid::ApplyOn::InnerIntersections<SGV>());
// contributions from periodic boundaries
- localizable_op.append(LocalAdvectionDgCouplingOperator<I, V, SGV, m, F, F, RGV, V>(*numerical_flux_),
+ localizable_op.append(LocalAdvectionDgCouplingOperator<I, V, SGV, m, F, F, RGV, V>(
+ local_mass_matrix_provider_, *numerical_flux_, /*compute_outside=*/false),
param,
- *(XT::Grid::ApplyOn::PeriodicBoundaryIntersectionsOnce<SGV>() && !(*periodicity_exception_)));
+ *(XT::Grid::ApplyOn::PeriodicBoundaryIntersections<SGV>() && !(*periodicity_exception_)));
// contributions from other boundaries by custom numerical flux
for (const auto& boundary_treatment : boundary_treatments_by_custom_numerical_flux_) {
const auto& boundary_op = *boundary_treatment.first;
@@ -155,53 +176,9 @@ protected:
const auto& filter = *boundary_treatment.second;
localizable_op.append(boundary_op, param, filter);
}
- } // ... append_standard_contributions(...)
-
- void
- append_local_mass_matrix_inversion(LocalizableOperatorBase<SGV, V, m, 1, F, SGV, m, 1, F, RGV, V>& localizable_op,
- RangeFunctionType& range_function) const
- {
- localizable_op.append(
- /*prepare=*/[]() {},
- /*apply_local=*/
- [&](const auto& element) {
- // (creating these objects before the grid walk and reusing them would be more efficient, but not thread safe)
- auto local_range = range_function.local_discrete_function(element);
- const auto& range_basis = local_range->basis();
- using E = XT::Grid::extract_entity_t<RGV>;
- const LocalElementIntegralBilinearForm<E, m, 1, F, F> local_l2_bilinear_form(
- LocalElementProductIntegrand<E, m, 1, F, F>(1.));
- local_range->dofs().assign_from(XT::LA::solve(
- XT::LA::convert_to<XT::LA::CommonDenseMatrix<F>>(local_l2_bilinear_form.apply2(range_basis, range_basis)),
- XT::LA::convert_to<XT::LA::CommonDenseVector<F>>(local_range->dofs()),
- {{"type", XT::LA::SolverOptions<XT::LA::CommonDenseMatrix<F>>::types().at(0)},
- {"post_check_solves_system", "-1"}}));
- },
- /*finalize=*/[]() {});
- } // ... append_local_mass_matrix_inversion(...)
-
-public:
- void apply(const VectorType& source, VectorType& range, const XT::Common::Parameter& param = {}) const override
- {
- // some checks
- DUNE_THROW_IF(!source.valid(), Exceptions::operator_error, "source contains inf or nan!");
- DUNE_THROW_IF(!(this->parameter_type() <= param.type()),
- Exceptions::operator_error,
- "this->parameter_type() = " << this->parameter_type() << "\n param.type() = " << param.type());
- range.set_all(0);
- auto source_function = make_discrete_function(this->source_space_, source);
- auto range_function = make_discrete_function(this->range_space_, range);
- // set up the actual operator
- auto localizable_op = make_localizable_operator(this->assembly_grid_view_, source_function, range_function);
- this->append_standard_contributions(localizable_op, param);
- // and apply it in a first grid walk
+ // and apply it in a grid walk
localizable_op.assemble(/*use_tbb=*/true);
DUNE_THROW_IF(!range.valid(), Exceptions::operator_error, "range contains inf or nan!");
- // prepare inversion of the mass matrix
- this->append_local_mass_matrix_inversion(localizable_op, range_function);
- // and actually do it in a second grid walk
- localizable_op.walk(/*use_tbb=*/true); // Do not call assemble() more than once, will not do anything!
- DUNE_THROW_IF(!range.valid(), Exceptions::operator_error, "range contains inf or nan!");
} // ... apply(...)
protected:
@@ -210,6 +187,7 @@ protected:
const SourceSpaceType& source_space_;
const RangeSpaceType& range_space_;
std::unique_ptr<XT::Grid::IntersectionFilter<SGV>> periodicity_exception_;
+ LocalMassMatrixProvider<RGV, m, 1, F> local_mass_matrix_provider_;
std::list<std::pair<std::unique_ptr<BoundaryTreatmentByCustomNumericalFluxOperatorType>,
std::unique_ptr<XT::Grid::IntersectionFilter<SGV>>>>
boundary_treatments_by_custom_numerical_flux_;