diff --git a/.gitignore b/.gitignore
index 197a93be4ec62814daac4395fd2872dfb58bd8b3..1701849adf8b3da45bffd6322cf43ef748e47749 100644
--- a/.gitignore
+++ b/.gitignore
@@ -46,3 +46,5 @@ demos-*
 build-*
 *.pyc
 .idea
+tags
+.vscode
diff --git a/dune/gdt/local/integrands/generic.hh b/dune/gdt/local/integrands/generic.hh
index ebf3744a14d06dfb3b39bce8a57cb899040e7b58..494170ebaf7ab2ab51aaea7ad2d10225cb4d358d 100644
--- a/dune/gdt/local/integrands/generic.hh
+++ b/dune/gdt/local/integrands/generic.hh
@@ -232,7 +232,7 @@ public:
     return std::make_unique<ThisType>(*this);
   }
 
-  virtual void post_bind(const I& intersect) const override final
+  virtual void post_bind(const I& intersect) override final
   {
     post_bind_(intersect);
   }
diff --git a/dune/gdt/local/operators/advection-dg.hh b/dune/gdt/local/operators/advection-dg.hh
index 3036e3a78dae8e76bf2b262db881261298f12320..1d316805dee6c002395a8281b8c58f7e4ad2a1b7 100644
--- a/dune/gdt/local/operators/advection-dg.hh
+++ b/dune/gdt/local/operators/advection-dg.hh
@@ -56,13 +56,27 @@ public:
   using FluxType = XT::Functions::FluxFunctionInterface<E, m, d, m, RF>;
   using LocalMassMatrixProviderType = LocalMassMatrixProvider<RGV, m, 1, RF>;
 
+  // When using this constructor, source has to be set by a call to with_source before calling apply
+  LocalAdvectionDgVolumeOperator(const FluxType& flux)
+    : BaseType(flux.parameter_type())
+    , flux_(flux)
+    , local_flux_(flux_.local_function())
+  {}
+
   LocalAdvectionDgVolumeOperator(const SourceType& source, const FluxType& flux)
     : BaseType(source, flux.parameter_type())
     , flux_(flux)
     , local_flux_(flux_.local_function())
   {}
 
-  /// Applies the inverse of the local mass matrix.
+  // When using this constructor, source has to be set by a call to with_source before calling apply
+  LocalAdvectionDgVolumeOperator(const LocalMassMatrixProviderType& local_mass_matrices, const FluxType& flux)
+    : BaseType(flux.parameter_type())
+    , flux_(flux)
+    , local_flux_(flux_.local_function())
+    , local_mass_matrices_(local_mass_matrices)
+  {}
+
   LocalAdvectionDgVolumeOperator(const SourceType& source,
                                  const LocalMassMatrixProviderType& local_mass_matrices,
                                  const FluxType& flux)
@@ -72,7 +86,6 @@ public:
     , local_mass_matrices_(local_mass_matrices)
   {}
 
-  /// Applies the inverse of the local mass matrix.
   LocalAdvectionDgVolumeOperator(const SourceSpaceType& source_space,
                                  const SV& source_vector,
                                  const LocalMassMatrixProviderType& local_mass_matrices,
@@ -178,6 +191,15 @@ public:
   using NumericalFluxType = NumericalFluxInterface<I, d, m, IRR>;
   using LocalMassMatrixProviderType = LocalMassMatrixProvider<IRGV, m, 1, IRR>;
 
+  // When using this constructor, source has to be set by a call to with_source before calling apply
+  LocalAdvectionDgCouplingOperator(const NumericalFluxType& numerical_flux, bool compute_outside = true)
+    : BaseType(numerical_flux.parameter_type())
+    , local_source_outside_(nullptr)
+    , numerical_flux_(numerical_flux.copy())
+    , local_flux_(numerical_flux_->flux().local_function())
+    , compute_outside_(compute_outside)
+  {}
+
   LocalAdvectionDgCouplingOperator(const SourceType& source,
                                    const NumericalFluxType& numerical_flux,
                                    bool compute_outside = true)
@@ -188,7 +210,18 @@ public:
     , compute_outside_(compute_outside)
   {}
 
-  /// Applies the inverse of the local mass matrix.
+  // When using this constructor, source has to be set by a call to with_source before calling apply
+  LocalAdvectionDgCouplingOperator(const LocalMassMatrixProviderType& local_mass_matrices,
+                                   const NumericalFluxType& numerical_flux,
+                                   bool compute_outside = true)
+    : BaseType(numerical_flux.parameter_type())
+    , local_source_outside_(nullptr)
+    , numerical_flux_(numerical_flux.copy())
+    , local_flux_(numerical_flux_->flux().local_function())
+    , compute_outside_(compute_outside)
+    , local_mass_matrices_(local_mass_matrices)
+  {}
+
   LocalAdvectionDgCouplingOperator(const SourceType& source,
                                    const LocalMassMatrixProviderType& local_mass_matrices,
                                    const NumericalFluxType& numerical_flux,
@@ -201,7 +234,6 @@ public:
     , local_mass_matrices_(local_mass_matrices)
   {}
 
-  /// Applies the inverse of the local mass matrix.
   LocalAdvectionDgCouplingOperator(const SourceSpaceType& source_space,
                                    const SV& source_vector,
                                    const LocalMassMatrixProviderType& local_mass_matrices,
@@ -358,7 +390,6 @@ public:
     , numerical_flux_order_(numerical_flux_order)
   {}
 
-  /// Applies the inverse of the local mass matrix.
   LocalAdvectionDgBoundaryTreatmentByCustomNumericalFluxOperator(
       const SourceType& source,
       const LocalMassMatrixProviderType& local_mass_matrices,
@@ -371,7 +402,6 @@ public:
     , local_mass_matrices_(local_mass_matrices)
   {}
 
-  /// Applies the inverse of the local mass matrix.
   LocalAdvectionDgBoundaryTreatmentByCustomNumericalFluxOperator(
       const SourceSpaceType& source_space,
       const SV& source_vector,
@@ -493,7 +523,6 @@ public:
     , extrapolate_(boundary_extrapolation_lambda)
   {}
 
-  /// Applies the inverse of the local mass matrix.
   LocalAdvectionDgBoundaryTreatmentByCustomExtrapolationOperator(
       const SourceType& source,
       const LocalMassMatrixProviderType& local_mass_matrices,
@@ -507,7 +536,6 @@ public:
     , local_mass_matrices_(local_mass_matrices)
   {}
 
-  /// Applies the inverse of the local mass matrix.
   LocalAdvectionDgBoundaryTreatmentByCustomExtrapolationOperator(
       const SourceSpaceType& source_space,
       const SV& source_vector,
@@ -608,6 +636,19 @@ public:
   using FluxType = XT::Functions::FunctionInterface<m, d, m, RF>;
   using LocalMassMatrixProviderType = LocalMassMatrixProvider<RGV, m, 1, RF>;
 
+  // When using this constructor, source has to be set by a call to with_source before calling apply
+  LocalAdvectionDgArtificialViscosityShockCapturingOperator(const SGV& assembly_grid_view,
+                                                            const double& nu_1 = 0.2,
+                                                            const double& alpha_1 = 1.0,
+                                                            const size_t index = 0)
+    : BaseType()
+    , local_source_outside_(nullptr)
+    , assembly_grid_view_(assembly_grid_view)
+    , nu_1_(nu_1)
+    , alpha_1_(alpha_1)
+    , index_(index)
+  {}
+
   LocalAdvectionDgArtificialViscosityShockCapturingOperator(const SourceType& source,
                                                             const SGV& assembly_grid_view,
                                                             const double& nu_1 = 0.2,
@@ -621,7 +662,21 @@ public:
     , index_(index)
   {}
 
-  /// Applies the inverse of the local mass matrix.
+  // When using this constructor, source has to be set by a call to with_source before calling apply
+  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()
+    , local_source_outside_(nullptr)
+    , assembly_grid_view_(assembly_grid_view)
+    , nu_1_(nu_1)
+    , alpha_1_(alpha_1)
+    , index_(index)
+    , local_mass_matrices_(local_mass_matrices)
+  {}
+
   LocalAdvectionDgArtificialViscosityShockCapturingOperator(const SourceType& source,
                                                             const LocalMassMatrixProviderType& local_mass_matrices,
                                                             const SGV& assembly_grid_view,
@@ -637,7 +692,6 @@ public:
     , local_mass_matrices_(local_mass_matrices)
   {}
 
-  /// Applies the inverse of the local mass matrix.
   LocalAdvectionDgArtificialViscosityShockCapturingOperator(
 
 
diff --git a/dune/gdt/local/operators/advection-fv.hh b/dune/gdt/local/operators/advection-fv.hh
index 3a2f3f00a2262771c21ed94aa339ff2d8b6bbd50..98202ba23ccecaba4453711d64d1e446130b79bd 100644
--- a/dune/gdt/local/operators/advection-fv.hh
+++ b/dune/gdt/local/operators/advection-fv.hh
@@ -64,6 +64,12 @@ public:
   using NumericalFluxType = NumericalFluxInterface<I, d, m, RR>;
   using LocalIntersectionCoords = typename NumericalFluxType::LocalIntersectionCoords;
 
+  // When using this constructor, source has to be set by a call to with_source before calling apply
+  LocalAdvectionFvCouplingOperator(const NumericalFluxType& numerical_flux, const bool source_x_independent = false)
+    : BaseType(numerical_flux.parameter_type())
+    , source_is_x_independent_(source_x_independent)
+  {}
+
   LocalAdvectionFvCouplingOperator(const SourceType& source,
                                    const NumericalFluxType& numerical_flux,
                                    const bool source_x_independent = false)
@@ -178,6 +184,16 @@ public:
   using LambdaType = std::function<StateType(
       const StateType& /*u*/, const StateDomainType& /*n*/, const XT::Common::Parameter& /*param*/)>;
 
+  // When using this constructor, source has to be set by a call to with_source before calling apply
+  LocalAdvectionFvBoundaryTreatmentByCustomNumericalFluxOperator(
+      LambdaType numerical_boundary_flux_lambda,
+      const XT::Common::ParameterType& boundary_treatment_param_type = {},
+      const bool source_is_x_independent = false)
+    : BaseType(boundary_treatment_param_type)
+    , numerical_boundary_flux_(numerical_boundary_flux_lambda)
+    , source_is_x_independent_(source_is_x_independent)
+  {}
+
   LocalAdvectionFvBoundaryTreatmentByCustomNumericalFluxOperator(
       const SourceType& source,
       LambdaType numerical_boundary_flux_lambda,
@@ -280,6 +296,18 @@ public:
                                              const StateType& /*u*/,
                                              const XT::Common::Parameter& /*param*/)>;
 
+  // When using this constructor, source has to be set by a call to with_source before calling apply
+  LocalAdvectionFvBoundaryTreatmentByCustomExtrapolationOperator(
+      const NumericalFluxType& numerical_flux,
+      LambdaType boundary_extrapolation_lambda,
+      const XT::Common::ParameterType& boundary_treatment_param_type = {},
+      const bool source_is_x_independent = false)
+    : BaseType(numerical_flux.parameter_type() + boundary_treatment_param_type)
+    , numerical_flux_(numerical_flux.copy())
+    , extrapolate_(boundary_extrapolation_lambda)
+    , source_is_x_independent_(source_is_x_independent)
+  {}
+
   LocalAdvectionFvBoundaryTreatmentByCustomExtrapolationOperator(
       const SourceType& source,
       const NumericalFluxType& numerical_flux,
diff --git a/dune/gdt/local/operators/interfaces.hh b/dune/gdt/local/operators/interfaces.hh
index 73b38cc65af5defd39f5c086c27890e0d49c731f..58a0a6c3aeb2a707bed659f5bd86786b62bfd7bc 100644
--- a/dune/gdt/local/operators/interfaces.hh
+++ b/dune/gdt/local/operators/interfaces.hh
@@ -71,6 +71,13 @@ public:
 
   using ThisType = LocalElementOperatorInterface<SV, SGV, s_r, s_rC, SR, r_r, r_rC, RR, RGV, RV>;
 
+  // Allows default construction, source has to be set by a call to with_source before calling apply
+  LocalElementOperatorInterface(const XT::Common::ParameterType& param_type = {})
+    : XT::Common::ParametricInterface(param_type)
+    , source_()
+    , local_source_(nullptr)
+  {}
+
   LocalElementOperatorInterface(const SourceType& source, const XT::Common::ParameterType& param_type = {})
     : XT::Common::ParametricInterface(param_type)
     , source_(source)
@@ -178,6 +185,13 @@ public:
   using ORGV = OutsideRangeGridView;
   using LocalOutsideRangeType = LocalDiscreteFunction<ORV, ORGV, r_r, r_rC, RF>;
 
+  // Allows default construction, source has to be set by a call to with_source before calling apply
+  LocalIntersectionOperatorInterface(const XT::Common::ParameterType& param_type = {})
+    : XT::Common::ParametricInterface(param_type)
+    , source_()
+    , local_source_(nullptr)
+  {}
+
   LocalIntersectionOperatorInterface(const SourceType& src, const XT::Common::ParameterType& param_type = {})
     : XT::Common::ParametricInterface(param_type)
     , source_(src)
diff --git a/dune/gdt/operators/interfaces.hh b/dune/gdt/operators/interfaces.hh
index bd0e31984dcc2a582541af10de21e9e08e537ab4..2737772388e1383226f8d602b4a4cb90792028eb 100644
--- a/dune/gdt/operators/interfaces.hh
+++ b/dune/gdt/operators/interfaces.hh
@@ -133,6 +133,7 @@ public:
   using F = FieldType;
 
   using SGV = SourceGridView;
+  using E = XT::Grid::extract_entity_t<SGV>;
   static const constexpr size_t s_r = source_dim;
   static const constexpr size_t s_rC = source_dim_cols;
 
@@ -141,6 +142,7 @@ public:
   static const constexpr size_t r_rC = range_dim_cols;
 
   using SourceSpaceType = SpaceInterface<SGV, s_r, s_rC, F>;
+  using SourceFunctionInterfaceType = XT::Functions::GridFunctionInterface<E, s_r, s_rC, F>;
   using SourceFunctionType = DiscreteFunction<V, SGV, s_r, s_rC, F>;
   using ConstSourceFunctionType = ConstDiscreteFunction<V, SGV, s_r, s_rC, F>;
 
diff --git a/dune/gdt/operators/localizable-operator.hh b/dune/gdt/operators/localizable-operator.hh
index 6fb63f9745e31db2ac304dac2f1c30eef922f0e2..d7d00e18eb259a238ce796a63a459fe7b3a3c803 100644
--- a/dune/gdt/operators/localizable-operator.hh
+++ b/dune/gdt/operators/localizable-operator.hh
@@ -222,6 +222,27 @@ make_localizable_operator_applicator(AGV assembly_grid_view,
       assembly_grid_view, source, range);
 }
 
+template <class AGV,
+          size_t s_r,
+          size_t s_rC,
+          class SF,
+          size_t r_r,
+          size_t r_rC,
+          class RF,
+          class RGV,
+          class RV,
+          class SV = RV>
+std::enable_if_t<XT::Grid::is_layer<AGV>::value,
+                 LocalizableDiscreteOperatorApplicator<AGV, SV, s_r, s_rC, SF, AGV, r_r, r_rC, RF, RGV, RV>>
+make_localizable_operator_applicator(
+    AGV assembly_grid_view,
+    const XT::Functions::GridFunctionInterface<XT::Grid::extract_entity_t<AGV>, s_r, s_rC, SF>& source,
+    DiscreteFunction<RV, RGV, r_r, r_rC, RF>& range)
+{
+  return LocalizableDiscreteOperatorApplicator<AGV, SV, s_r, s_rC, SF, AGV, r_r, r_rC, RF, RGV, RV>(
+      assembly_grid_view, source, range);
+}
+
 template <class AGV,
           class SV,
           size_t s_r,
@@ -273,6 +294,7 @@ public:
 
   using typename BaseType::MatrixOperatorType;
   using typename BaseType::RangeSpaceType;
+  using typename BaseType::SourceFunctionInterfaceType;
   using typename BaseType::SourceSpaceType;
   using typename BaseType::VectorType;
 
@@ -326,36 +348,42 @@ public:
 
   using BaseType::apply;
 
-  void apply(const VectorType& source, VectorType& range, const XT::Common::Parameter& param = {}) const override
+  void apply(const SourceFunctionInterfaceType& source_function,
+             VectorType& range,
+             const XT::Common::Parameter& param = {}) const
   {
-    // 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_applicator(this->assembly_grid_view_, source_function, range_function);
     // - element contributions
     for (const auto& op_and_filter : local_element_operators_) {
-      const auto& local_op = *op_and_filter.first;
+      const auto local_op = op_and_filter.first->with_source(source_function);
       const auto& filter = *op_and_filter.second;
-      localizable_op.append(local_op, param, filter);
+      localizable_op.append(*local_op, param, filter);
     }
     // - intersection contributions
     for (const auto& op_and_filter : local_intersection_operators_) {
-      const auto& local_op = *op_and_filter.first;
+      const auto local_op = op_and_filter.first->with_source(source_function);
       const auto& filter = *op_and_filter.second;
-      localizable_op.append(local_op, param, filter);
+      localizable_op.append(*local_op, param, filter);
     }
     // 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!");
   } // ... apply(...)
 
+  void apply(const VectorType& source, VectorType& range, const XT::Common::Parameter& param = {}) const override
+  {
+    DUNE_THROW_IF(!source.valid(), Exceptions::operator_error, "source contains inf or nan!");
+    const auto source_function = make_discrete_function(this->source_space_, source);
+    apply(source_function, range, param);
+  } // ... apply(...)
+
   std::vector<std::string> jacobian_options() const override final
   {
     return {"finite-differences"};
@@ -386,16 +414,17 @@ public:
     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
+    const auto source_function = make_discrete_function(this->source_space_, source);
     for (const auto& op_and_filter : local_element_operators_) {
-      const auto& local_op = *op_and_filter.first;
+      const auto local_op = op_and_filter.first->with_source(source_function);
       const auto& filter = *op_and_filter.second;
-      jacobian_op.append(local_op, source, parameter, filter);
+      jacobian_op.append(*local_op, source, parameter, filter);
     }
     // - intersection contributions
     for (const auto& op_and_filter : local_intersection_operators_) {
-      const auto& local_op = *op_and_filter.first;
+      const auto local_op = op_and_filter.first->with_source(source_function);
       const auto& filter = *op_and_filter.second;
-      jacobian_op.append(local_op, source, parameter, filter);
+      jacobian_op.append(*local_op, source, parameter, filter);
     }
   } // ... jacobian(...)