Skip to content
Snippets Groups Projects
Commit 6cddf6ac authored by Dr. Felix Tobias Schindler's avatar Dr. Felix Tobias Schindler
Browse files

[operators.advection-dg] add artificial viscosity, jacobian

parent 04f605f2
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/gdt/local/operators/advection-dg.hh
+ 143
− 0
View file @ 6cddf6ac
...@@ -14,9 +14,11 @@ ...@@ -14,9 +14,11 @@
#include <functional> #include <functional>
#include <dune/geometry/quadraturerules.hh> #include <dune/geometry/quadraturerules.hh>
#include <dune/grid/common/rangegenerators.hh>
#include <dune/xt/common/memory.hh> #include <dune/xt/common/memory.hh>
#include <dune/xt/common/parameter.hh> #include <dune/xt/common/parameter.hh>
#include <dune/xt/grid/intersection.hh>
#include <dune/gdt/exceptions.hh> #include <dune/gdt/exceptions.hh>
#include <dune/gdt/local/dof-vector.hh> #include <dune/gdt/local/dof-vector.hh>
...@@ -461,6 +463,147 @@ private: ...@@ -461,6 +463,147 @@ private:
}; // class LocalAdvectionDgBoundaryTreatmentByCustomExtrapolationOperator }; // class LocalAdvectionDgBoundaryTreatmentByCustomExtrapolationOperator
template <class SV, class SGV, size_t m = 1, class SF = double, class RF = SF, class RGV = SGV, class RV = SV>
class LocalAdvectionDgArtificialViscosityShockCapturingOperator
: public LocalElementOperatorInterface<SV, SGV, m, 1, SF, m, 1, RF, RGV, RV>
{
using ThisType = LocalAdvectionDgArtificialViscosityShockCapturingOperator<SV, SGV, m, SF, RF, RGV, RV>;
using BaseType = LocalElementOperatorInterface<SV, SGV, m, 1, SF, m, 1, RF, RGV, RV>;
public:
using BaseType::d;
using typename BaseType::D;
using typename BaseType::LocalRangeType;
using typename BaseType::SourceType;
using FluxType = XT::Functions::FunctionInterface<m, d, m, RF>;
using LocalMassMatrixProviderType = LocalMassMatrixProvider<RGV, m, 1, RF>;
LocalAdvectionDgArtificialViscosityShockCapturingOperator(const SGV& assembly_grid_view,
const double& nu_1 = 0.2,
const double& alpha_1 = 1.0,
const size_t index = 0)
: BaseType()
, assembly_grid_view_(assembly_grid_view)
, nu_1_(nu_1)
, alpha_1_(alpha_1)
, index_(index)
{}
/// Applies the inverse of the local mass matrix.
LocalAdvectionDgArtificialViscosityShockCapturingOperator(const LocalMassMatrixProviderType& local_mass_matrices,
const SGV& assembly_grid_view,
const double& nu_1 = 0.2,
const double& alpha_1 = 1.0,
const size_t index = 0)
: BaseType()
, assembly_grid_view_(assembly_grid_view)
, nu_1_(nu_1)
, alpha_1_(alpha_1)
, index_(index)
, local_mass_matrices_(local_mass_matrices)
{}
LocalAdvectionDgArtificialViscosityShockCapturingOperator(const ThisType& other)
: BaseType(other)
, assembly_grid_view_(other.assembly_grid_view_)
, nu_1_(other.nu_1_)
, alpha_1_(other.alpha_1_)
, index_(other.index_)
, local_mass_matrices_(other.local_mass_matrices_)
{}
std::unique_ptr<BaseType> copy() const override final
{
return std::make_unique<ThisType>(*this);
}
void apply(const SourceType& source,
LocalRangeType& local_range,
const XT::Common::Parameter& param = {}) const override final
{
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_element = source.local_discrete_function(element);
if (local_source_element->order(param) <= 0 || basis.order(param) <= 0)
return;
// compute jump indicator (8.176)
double element_jump_indicator = 0;
double element_boundary_without_domain_boundary = (d == 1) ? 1. : 0.;
for (auto&& intersection : intersections(assembly_grid_view_, element)) {
if (intersection.neighbor() && !intersection.boundary()) {
if (d > 1)
element_boundary_without_domain_boundary += XT::Grid::diameter(intersection);
const auto neighbor = intersection.outside();
const auto local_source_neighbor = source.local_discrete_function(neighbor);
const auto integration_order =
std::pow(std::max(local_source_element->order(param), local_source_neighbor->order(param)), 2);
for (auto&& quadrature_point :
QuadratureRules<D, d - 1>::rule(intersection.geometry().type(), integration_order)) {
const auto point_in_reference_intersection = quadrature_point.position();
const auto point_in_reference_element =
intersection.geometryInInside().global(point_in_reference_intersection);
const auto point_in_reference_neighbor =
intersection.geometryInOutside().global(point_in_reference_intersection);
const auto integration_factor = intersection.geometry().integrationElement(point_in_reference_intersection);
const auto quadrature_weight = quadrature_point.weight();
const auto value_on_element = local_source_element->evaluate(point_in_reference_element, param)[index_];
const auto value_on_neighbor = local_source_neighbor->evaluate(point_in_reference_neighbor, param)[index_];
element_jump_indicator +=
integration_factor * quadrature_weight * std::pow(value_on_element - value_on_neighbor, 2);
}
}
element_jump_indicator /= element_boundary_without_domain_boundary * element.geometry().volume();
}
// compute smoothed discrete jump indicator (8.180)
double smoothed_discrete_jump_indicator = 0;
const double xi_min = 0.5;
const double xi_max = 1.5;
if (element_jump_indicator < xi_min)
smoothed_discrete_jump_indicator = 0;
else if (!(element_jump_indicator < xi_max))
smoothed_discrete_jump_indicator = 1;
else
smoothed_discrete_jump_indicator =
0.5 * std::sin(M_PI * (element_jump_indicator - (xi_max - xi_min)) / (2 * (xi_max - xi_min))) + 0.5;
// evaluate artificial viscosity form (8.183)
if (smoothed_discrete_jump_indicator > 0) {
const auto h = element.geometry().volume();
for (const auto& quadrature_point : QuadratureRules<D, d>::rule(
element.type(),
std::max(0, local_source_element->order(param) - 1) + std::max(0, basis.order(param) - 1))) {
const auto point_in_reference_element = quadrature_point.position();
const auto integration_factor = element.geometry().integrationElement(point_in_reference_element);
const auto quadrature_weight = quadrature_point.weight();
const auto source_jacobian = local_source_element->jacobian(point_in_reference_element, param);
basis.jacobians(point_in_reference_element, basis_jacobians_, param);
// compute beta_h
for (size_t ii = 0; ii < basis.size(param); ++ii)
local_dofs_[ii] += integration_factor * quadrature_weight * nu_1_ * std::pow(h, alpha_1_)
* smoothed_discrete_jump_indicator * (source_jacobian[0] * basis_jacobians_[ii][0]);
}
}
// 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 SGV& assembly_grid_view_;
const double nu_1_;
const double alpha_1_;
const size_t index_;
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 LocalAdvectionDgArtificialViscosityShockCapturingOperator
} // namespace GDT } // namespace GDT
} // namespace Dune } // namespace Dune
......
dune/gdt/operators/advection-dg.hh
+ 104
− 171
View file @ 6cddf6ac
...@@ -38,6 +38,22 @@ namespace Dune { ...@@ -38,6 +38,22 @@ namespace Dune {
namespace GDT { namespace GDT {
static double advection_dg_artificial_viscosity_default_nu_1()
{
return 0.2;
}
static double advection_dg_artificial_viscosity_default_alpha_1()
{
return 1.0;
}
static size_t advection_dg_artificial_viscosity_default_component()
{
return 0;
}
/** /**
* \note See OperatorInterface for a description of the template arguments. * \note See OperatorInterface for a description of the template arguments.
* *
...@@ -65,6 +81,7 @@ public: ...@@ -65,6 +81,7 @@ public:
using BoundaryTreatmentByCustomExtrapolationOperatorType = using BoundaryTreatmentByCustomExtrapolationOperatorType =
LocalAdvectionDgBoundaryTreatmentByCustomExtrapolationOperator<I, V, SGV, m, F, F, RGV, V>; LocalAdvectionDgBoundaryTreatmentByCustomExtrapolationOperator<I, V, SGV, m, F, F, RGV, V>;
using typename BaseType::MatrixOperatorType;
using typename BaseType::RangeFunctionType; using typename BaseType::RangeFunctionType;
using typename BaseType::RangeSpaceType; using typename BaseType::RangeSpaceType;
using typename BaseType::SourceSpaceType; using typename BaseType::SourceSpaceType;
...@@ -75,7 +92,10 @@ public: ...@@ -75,7 +92,10 @@ public:
const NumericalFluxType& numerical_flux, const NumericalFluxType& numerical_flux,
const SourceSpaceType& source_space, const SourceSpaceType& source_space,
const RangeSpaceType& range_space, const RangeSpaceType& range_space,
const XT::Grid::IntersectionFilter<SGV>& periodicity_exception = XT::Grid::ApplyOn::NoIntersections<SGV>()) const XT::Grid::IntersectionFilter<SGV>& periodicity_exception = XT::Grid::ApplyOn::NoIntersections<SGV>(),
const double& artificial_viscosity_nu_1 = advection_dg_artificial_viscosity_default_nu_1(),
const double& artificial_viscosity_alpha_1 = advection_dg_artificial_viscosity_default_alpha_1(),
const size_t artificial_viscosity_component = advection_dg_artificial_viscosity_default_component())
: BaseType(numerical_flux.parameter_type()) : BaseType(numerical_flux.parameter_type())
, assembly_grid_view_(assembly_grid_view) , assembly_grid_view_(assembly_grid_view)
, numerical_flux_(numerical_flux.copy()) , numerical_flux_(numerical_flux.copy())
...@@ -83,6 +103,9 @@ public: ...@@ -83,6 +103,9 @@ public:
, range_space_(range_space) , range_space_(range_space)
, periodicity_exception_(periodicity_exception.copy()) , periodicity_exception_(periodicity_exception.copy())
, local_mass_matrix_provider_(assembly_grid_view_, range_space_) , local_mass_matrix_provider_(assembly_grid_view_, range_space_)
, artificial_viscosity_nu_1_(artificial_viscosity_nu_1)
, artificial_viscosity_alpha_1_(artificial_viscosity_alpha_1)
, artificial_viscosity_component_(artificial_viscosity_component)
{ {
// we assemble these once, to be used in each apply later on // we assemble these once, to be used in each apply later on
auto walker = XT::Grid::make_walker(assembly_grid_view_); auto walker = XT::Grid::make_walker(assembly_grid_view_);
...@@ -176,11 +199,88 @@ public: ...@@ -176,11 +199,88 @@ public:
const auto& filter = *boundary_treatment.second; const auto& filter = *boundary_treatment.second;
localizable_op.append(boundary_op, param, filter); localizable_op.append(boundary_op, param, filter);
} }
// artificial viscosity by shock capturing [DF2015, Sec. 8.5]
localizable_op.append(LocalAdvectionDgArtificialViscosityShockCapturingOperator<V, SGV, m, F, F, RGV, V>(
local_mass_matrix_provider_,
this->assembly_grid_view_,
artificial_viscosity_nu_1_,
artificial_viscosity_alpha_1_,
artificial_viscosity_component_),
param);
// and apply it in a grid walk // and apply it in a grid walk
localizable_op.assemble(/*use_tbb=*/true); localizable_op.assemble(/*use_tbb=*/true);
DUNE_THROW_IF(!range.valid(), Exceptions::operator_error, "range contains inf or nan!"); DUNE_THROW_IF(!range.valid(), Exceptions::operator_error, "range contains inf or nan!");
} // ... apply(...) } // ... apply(...)
std::vector<std::string> jacobian_options() const override final
{
return {"finite-differences"};
}
XT::Common::Configuration jacobian_options(const std::string& type) const override final
{
DUNE_THROW_IF(type != this->jacobian_options().at(0), Exceptions::operator_error, "type = " << type);
return {{"type", type}, {"eps", "1e-7"}};
}
using BaseType::jacobian;
void jacobian(const VectorType& source,
MatrixOperatorType& jacobian_op,
const XT::Common::Configuration& opts,
const XT::Common::Parameter& param = {}) const override final
{
// 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());
DUNE_THROW_IF(!opts.has_key("type"), Exceptions::operator_error, opts);
DUNE_THROW_IF(opts.get<std::string>("type") != jacobian_options().at(0), Exceptions::operator_error, opts);
const auto default_opts = jacobian_options(jacobian_options().at(0));
const auto eps = opts.get("eps", default_opts.template get<double>("eps"));
const auto parameter = param + XT::Common::Parameter({"finite-difference-jacobians.eps", eps});
// append the same local ops with the same filters as in apply() above
// element contributions
jacobian_op.append(
LocalAdvectionDgVolumeOperator<V, SGV, m, F, F, RGV, V>(local_mass_matrix_provider_, numerical_flux_->flux()),
source,
parameter);
// contributions from inner intersections
jacobian_op.append(LocalAdvectionDgCouplingOperator<I, V, SGV, m, F, F, RGV, V>(
local_mass_matrix_provider_, *numerical_flux_, /*compute_outside=*/false),
source,
parameter,
XT::Grid::ApplyOn::InnerIntersections<SGV>());
// contributions from periodic boundaries
jacobian_op.append(LocalAdvectionDgCouplingOperator<I, V, SGV, m, F, F, RGV, V>(
local_mass_matrix_provider_, *numerical_flux_, /*compute_outside=*/false),
source,
parameter,
*(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;
const auto& filter = *boundary_treatment.second;
jacobian_op.append(boundary_op, source, parameter, filter);
}
// contributions from other boundaries by custom extrapolation
for (const auto& boundary_treatment : boundary_treatments_by_custom_extrapolation_) {
const auto& boundary_op = *boundary_treatment.first;
const auto& filter = *boundary_treatment.second;
jacobian_op.append(boundary_op, source, parameter, filter);
}
// artificial viscosity by shock capturing [DF2015, Sec. 8.5]
jacobian_op.append(LocalAdvectionDgArtificialViscosityShockCapturingOperator<V, SGV, m, F, F, RGV, V>(
local_mass_matrix_provider_,
this->assembly_grid_view_,
artificial_viscosity_nu_1_,
artificial_viscosity_alpha_1_,
artificial_viscosity_component_),
source,
param);
} // ... jacobian(...)
protected: protected:
const SGV assembly_grid_view_; const SGV assembly_grid_view_;
const std::unique_ptr<const NumericalFluxType> numerical_flux_; const std::unique_ptr<const NumericalFluxType> numerical_flux_;
...@@ -188,6 +288,9 @@ protected: ...@@ -188,6 +288,9 @@ protected:
const RangeSpaceType& range_space_; const RangeSpaceType& range_space_;
std::unique_ptr<XT::Grid::IntersectionFilter<SGV>> periodicity_exception_; std::unique_ptr<XT::Grid::IntersectionFilter<SGV>> periodicity_exception_;
LocalMassMatrixProvider<RGV, m, 1, F> local_mass_matrix_provider_; LocalMassMatrixProvider<RGV, m, 1, F> local_mass_matrix_provider_;
const double artificial_viscosity_nu_1_;
const double artificial_viscosity_alpha_1_;
const size_t artificial_viscosity_component_;
std::list<std::pair<std::unique_ptr<BoundaryTreatmentByCustomNumericalFluxOperatorType>, std::list<std::pair<std::unique_ptr<BoundaryTreatmentByCustomNumericalFluxOperatorType>,
std::unique_ptr<XT::Grid::IntersectionFilter<SGV>>>> std::unique_ptr<XT::Grid::IntersectionFilter<SGV>>>>
boundary_treatments_by_custom_numerical_flux_; boundary_treatments_by_custom_numerical_flux_;
...@@ -211,176 +314,6 @@ make_advection_dg_operator( ...@@ -211,176 +314,6 @@ make_advection_dg_operator(
} }
#if 0
/**
* \note See AdvectionDgOperator for a description of the template arguments.
*
* \sa AdvectionDgOperator
* \sa [FK2007, Sec. 6.2] for nu_1
*/
template <class M, class SGV, size_t m = 1, class RGV = SGV>
class AdvectionDgArtificialViscosityOperator : public AdvectionDgOperator<M, SGV, m, RGV>
{
using ThisType = AdvectionDgArtificialViscosityOperator<M, SGV, m, RGV>;
using BaseType = AdvectionDgOperator<M, SGV, m, RGV>;
using typename BaseType::D;
static constexpr size_t d = BaseType::d;
public:
using typename BaseType::F;
using typename BaseType::NumericalFluxType;
using typename BaseType::RangeSpaceType;
using typename BaseType::SourceSpaceType;
using typename BaseType::VectorType;
AdvectionDgArtificialViscosityOperator(
const SGV& assembly_grid_view,
const NumericalFluxType& numerical_flux,
const SourceSpaceType& source_space,
const RangeSpaceType& range_space,
const XT::Grid::IntersectionFilter<SGV>& periodicity_exception = XT::Grid::ApplyOn::NoIntersections<SGV>(),
const double nu_1 = 0.2,
const double alpha_1 = 1.,
const size_t jump_indicator_component = 0)
: BaseType(assembly_grid_view, numerical_flux, source_space, range_space, periodicity_exception)
, jump_indicator_component_(jump_indicator_component)
, nu_1_(nu_1)
, alpha_1_(alpha_1)
{}
AdvectionDgArtificialViscosityOperator(ThisType&& source) = default;
using BaseType::apply;
void apply(const VectorType& source, VectorType& range, const XT::Common::Parameter& param = {}) const override final
{
// 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);
# if 0
auto src = source.copy();
auto source_function = make_discrete_function(source_space_, src);
// apply local P0 projection
auto walker = XT::Grid::make_walker(assembly_grid_view_);
walker.append([&](const auto& element) {
auto local_source = source_function.local_discrete_function(element);
using E = XT::Grid::extract_entity_t<SGV>;
const auto average = LocalElementIntegralFunctional<E>(LocalElementIdentityIntegrand<E>()).apply(*local_source)[0]
/ element.geometry().volume();
local_source->dofs().assign_from(source_space_.finite_element(element.geometry().type())
.interpolation()
.interpolate([&](const auto& /*xx*/) { return average; }, 0));
});
walker.walk(/*use_tbb=*/true);
# endif
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);
// compute jump indicators for artificial viscosity at detected shocks, see DF2015, p. 449, (8.180)
// we need a thread safe vector with one entry per grid element and just use a scalar FV discrete function
const auto fv_space = make_finite_volume_space<1>(this->assembly_grid_view_);
auto jump_indicators = make_discrete_function<XT::LA::CommonDenseVector<F>>(fv_space);
localizable_op.append(
/*prepare=*/[]() {},
/*apply_local=*/
[&](const auto& element) {
// compute jump indicator (8.176)
double element_jump_indicator = 0;
double element_boundary_without_domain_boundary = (d == 1) ? 1. : 0.;
const auto local_source_element = source_function.local_discrete_function(element);
for (auto&& intersection : intersections(this->assembly_grid_view_, element)) {
if (intersection.neighbor() && !intersection.boundary()) {
if (d > 1)
element_boundary_without_domain_boundary += XT::Grid::diameter(intersection);
const auto neighbor = intersection.outside();
const auto local_source_neighbor = source_function.local_discrete_function(neighbor);
const auto integration_order =
std::pow(std::max(local_source_element->order(), local_source_neighbor->order()), 2);
for (auto&& quadrature_point :
QuadratureRules<D, d - 1>::rule(intersection.geometry().type(), integration_order)) {
const auto point_in_reference_intersection = quadrature_point.position();
const auto point_in_reference_element =
intersection.geometryInInside().global(point_in_reference_intersection);
const auto point_in_reference_neighbor =
intersection.geometryInOutside().global(point_in_reference_intersection);
const auto integration_factor =
intersection.geometry().integrationElement(point_in_reference_intersection);
const auto quadrature_weight = quadrature_point.weight();
const auto value_on_element =
local_source_element->evaluate(point_in_reference_element)[jump_indicator_component_];
const auto value_on_neighbor =
local_source_neighbor->evaluate(point_in_reference_neighbor)[jump_indicator_component_];
element_jump_indicator +=
integration_factor * quadrature_weight * std::pow(value_on_element - value_on_neighbor, 2);
}
}
element_jump_indicator /= element_boundary_without_domain_boundary * element.geometry().volume();
}
// compute smoothed discrete jump indicator (8.180)
double smoothed_discrete_jump_indicator = 0;
const double xi_min = 0.5;
const double xi_max = 1.5;
if (element_jump_indicator < xi_min)
smoothed_discrete_jump_indicator = 0;
else if (!(element_jump_indicator < xi_max))
smoothed_discrete_jump_indicator = 1;
else
smoothed_discrete_jump_indicator =
0.5 * std::sin(M_PI * (element_jump_indicator - (xi_max - xi_min)) / (2 * (xi_max - xi_min))) + 0.5;
jump_indicators.local_discrete_function(element)->dofs()[0] = smoothed_discrete_jump_indicator;
},
/*finalize=*/[]() {});
// do the actual (first) grid walk: the above operators will be applied and afterwards cleared
localizable_op.assemble(/*use_tbb=*/true);
DUNE_THROW_IF(!range.valid(), Exceptions::operator_error, "range contains inf or nan!");
// compute artificial viscosity shock capturing [see DF2015, p. 450, (8.183, 8.1.84)] by appending a grid element
// functor, use the localizable_op as a XT::Grid::Walker
localizable_op.append(
/*prepare=*/[]() {},
/*apply_local=*/
[&](const auto& element) {
// evaluate artificial viscosity form (8.183)
const auto local_source = source_function.local_discrete_function(element);
auto local_range = range_function.local_discrete_function(element);
const auto& local_basis = local_range->basis();
const auto h = element.geometry().volume();
const auto jump_indicator_element = jump_indicators.local_discrete_function(element)->dofs()[0];
for (const auto& quadrature_point : QuadratureRules<D, d>::rule(element.type(), 2 * local_basis.order())) {
const auto point_in_reference_element = quadrature_point.position();
const auto integration_factor = element.geometry().integrationElement(point_in_reference_element);
const auto quadrature_weight = quadrature_point.weight();
const auto source_jacobian = local_source->jacobian(point_in_reference_element);
const auto basis_jacobians = local_basis.jacobians_of_set(point_in_reference_element);
// compute beta_h
for (size_t ii = 0; ii < local_basis.size(); ++ii)
local_range->dofs()[ii] += integration_factor * quadrature_weight * nu_1_ * std::pow(h, alpha_1_)
* jump_indicator_element * (source_jacobian[0] * basis_jacobians[ii][0]);
}
},
/*finalize=*/[]() {});
// do the actual (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 the inverse mass matrix by appending a grid element functor, use the localizable_op as a XT::Grid::Walker
this->append_local_mass_matrix_inversion(localizable_op, range_function);
// do the actual (third) 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(...)
private:
const size_t jump_indicator_component_;
const double nu_1_;
const double alpha_1_;
}; // class AdvectionDgArtificialViscosityOperator
#endif
} // namespace GDT } // namespace GDT
} // namespace Dune } // namespace Dune
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment