diff --git a/dumux/assembly/boxlocalresidual.hh b/dumux/assembly/boxlocalresidual.hh
index 72eec7960c8b80853b4f9a9fb2e6b53cbc648825..cb8cdf25ab9f54875aa28d559c74c70c753777a5 100644
--- a/dumux/assembly/boxlocalresidual.hh
+++ b/dumux/assembly/boxlocalresidual.hh
@@ -52,7 +52,7 @@ class BoxLocalResidual : public FVLocalResidual<TypeTag>
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
     using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
 
 public:
     using ElementResidualVector = typename ParentType::ElementResidualVector;
@@ -84,15 +84,15 @@ public:
     }
 
     //! evaluate flux residuals for one sub control volume face
-    ResidualVector evalFlux(const Problem& problem,
-                            const Element& element,
-                            const FVElementGeometry& fvGeometry,
-                            const ElementVolumeVariables& elemVolVars,
-                            const ElementBoundaryTypes& elemBcTypes,
-                            const ElementFluxVariablesCache& elemFluxVarsCache,
-                            const SubControlVolumeFace& scvf) const
+    NumEqVector evalFlux(const Problem& problem,
+                         const Element& element,
+                         const FVElementGeometry& fvGeometry,
+                         const ElementVolumeVariables& elemVolVars,
+                         const ElementBoundaryTypes& elemBcTypes,
+                         const ElementFluxVariablesCache& elemFluxVarsCache,
+                         const SubControlVolumeFace& scvf) const
     {
-        ResidualVector flux(0.0);
+        NumEqVector flux(0.0);
 
         // inner faces
         if (!scvf.boundary())
diff --git a/dumux/assembly/cclocalassembler.hh b/dumux/assembly/cclocalassembler.hh
index 314b5702c7363983074f241108193907dd79f1b8..7e3f2ec3f99877923f20e024d33e52bbd5554c93 100644
--- a/dumux/assembly/cclocalassembler.hh
+++ b/dumux/assembly/cclocalassembler.hh
@@ -61,7 +61,7 @@ class CCLocalAssemblerBase : public FVLocalAssemblerBase<TypeTag, Assembler, Imp
     using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
-    using LocalResidualValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
 
 public:
 
@@ -133,7 +133,7 @@ class CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/tru
     using ThisType = CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, true>;
     using ParentType = CCLocalAssemblerBase<TypeTag, Assembler, ThisType, true>;
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-    using LocalResidualValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using FVElementGeometry = typename FVGridGeometry::LocalView;
@@ -157,7 +157,7 @@ public:
      *
      * \return The element residual at the current solution.
      */
-    LocalResidualValues assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
+    NumEqVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
     {
         //////////////////////////////////////////////////////////////////////////////////////////////////
         // Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the //
@@ -182,7 +182,7 @@ public:
         neighborElements.resize(numNeighbors);
 
         // assemble the undeflected residual
-        using Residuals = ReservedBlockVector<LocalResidualValues, maxElementStencilSize>;
+        using Residuals = ReservedBlockVector<NumEqVector, maxElementStencilSize>;
         Residuals origResiduals(numNeighbors + 1); origResiduals = 0.0;
         origResiduals[0] = this->evalLocalResidual()[0];
 
@@ -311,7 +311,7 @@ class CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/fal
     using ThisType = CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, false>;
     using ParentType = CCLocalAssemblerBase<TypeTag, Assembler, ThisType, false>;
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-    using LocalResidualValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
     using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
     using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
@@ -327,7 +327,7 @@ public:
      *
      * \return The element residual at the current solution.
      */
-    LocalResidualValues assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
+    NumEqVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
     {
         if (this->assembler().isStationaryProblem())
             DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
@@ -361,7 +361,7 @@ public:
         // element solution container to be deflected
         auto elemSol = elementSolution(element, curSol, fvGridGeometry);
 
-        LocalResidualValues partialDeriv;
+        NumEqVector partialDeriv;
 
         // derivatives in the neighbors with repect to the current elements
         for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
@@ -421,7 +421,7 @@ class CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/tr
 {
     using ThisType = CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, true>;
     using ParentType = CCLocalAssemblerBase<TypeTag, Assembler, ThisType, true>;
-    using LocalResidualValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
     using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
 
@@ -434,7 +434,7 @@ public:
      *
      * \return The element residual at the current solution.
      */
-    LocalResidualValues assembleJacobianAndResidualImpl(JacobianMatrix& A, const GridVariables& gridVariables)
+    NumEqVector assembleJacobianAndResidualImpl(JacobianMatrix& A, const GridVariables& gridVariables)
     {
         // assemble the undeflected residual
         const auto residual = this->evalLocalResidual()[0];
@@ -500,7 +500,7 @@ class CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/fa
 {
     using ThisType = CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, false>;
     using ParentType = CCLocalAssemblerBase<TypeTag, Assembler, ThisType, false>;
-    using LocalResidualValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
     using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
 
@@ -513,7 +513,7 @@ public:
      *
      * \return The element residual at the current solution.
      */
-    LocalResidualValues assembleJacobianAndResidualImpl(JacobianMatrix& A, const GridVariables& gridVariables)
+    NumEqVector assembleJacobianAndResidualImpl(JacobianMatrix& A, const GridVariables& gridVariables)
     {
         // assemble the undeflected residual
         const auto residual = this->evalLocalResidual()[0];
diff --git a/dumux/assembly/cclocalresidual.hh b/dumux/assembly/cclocalresidual.hh
index a94b4bc9ae6040175447b6443a4ddb4abcbec429..f8e5995097de9349941dc589699e5cea763f1cec 100644
--- a/dumux/assembly/cclocalresidual.hh
+++ b/dumux/assembly/cclocalresidual.hh
@@ -42,7 +42,7 @@ class CCLocalResidual : public FVLocalResidual<TypeTag>
     using ParentType = FVLocalResidual<TypeTag>;
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
@@ -69,14 +69,14 @@ public:
     }
 
     //! evaluate the flux residual for a sub control volume face
-    ResidualVector evalFlux(const Problem& problem,
-                            const Element& element,
-                            const FVElementGeometry& fvGeometry,
-                            const ElementVolumeVariables& elemVolVars,
-                            const ElementFluxVariablesCache& elemFluxVarsCache,
-                            const SubControlVolumeFace& scvf) const
+    NumEqVector evalFlux(const Problem& problem,
+                         const Element& element,
+                         const FVElementGeometry& fvGeometry,
+                         const ElementVolumeVariables& elemVolVars,
+                         const ElementFluxVariablesCache& elemFluxVarsCache,
+                         const SubControlVolumeFace& scvf) const
     {
-        ResidualVector flux(0.0);
+        NumEqVector flux(0.0);
 
         // inner faces
         if (!scvf.boundary())
diff --git a/dumux/assembly/fvlocalresidual.hh b/dumux/assembly/fvlocalresidual.hh
index 3a081c9ff1f2fcafd11af4180aa3b62efcdee606..5e1ecde657f27299d1bc66611ee51d8b07d3de41 100644
--- a/dumux/assembly/fvlocalresidual.hh
+++ b/dumux/assembly/fvlocalresidual.hh
@@ -54,7 +54,7 @@ class FVLocalResidual
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
     using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
     using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
@@ -64,7 +64,7 @@ class FVLocalResidual
 
 public:
     //! the container storing all element residuals
-    using ElementResidualVector = ReservedBlockVector<ResidualVector, FVElementGeometry::maxNumElementScvs>;
+    using ElementResidualVector = ReservedBlockVector<NumEqVector, FVElementGeometry::maxNumElementScvs>;
 
     //! the constructor
     FVLocalResidual(const Problem* problem,
@@ -269,9 +269,9 @@ public:
      * \note has to be implemented by the model specific residual class
      *
      */
-    ResidualVector computeStorage(const Problem& problem,
-                                  const SubControlVolume& scv,
-                                  const VolumeVariables& volVars) const
+    NumEqVector computeStorage(const Problem& problem,
+                               const SubControlVolume& scv,
+                               const VolumeVariables& volVars) const
     {
         DUNE_THROW(Dune::NotImplemented, "This model does not implement a storage method!");
     }
@@ -289,13 +289,13 @@ public:
      *       in the user interface of the problem
      *
      */
-    ResidualVector computeSource(const Problem& problem,
-                                 const Element& element,
-                                 const FVElementGeometry& fvGeometry,
-                                 const ElementVolumeVariables& elemVolVars,
-                                 const SubControlVolume &scv) const
+    NumEqVector computeSource(const Problem& problem,
+                              const Element& element,
+                              const FVElementGeometry& fvGeometry,
+                              const ElementVolumeVariables& elemVolVars,
+                              const SubControlVolume &scv) const
     {
-        ResidualVector source(0.0);
+        NumEqVector source(0.0);
 
         // add contributions from volume flux sources
         source += problem.source(element, fvGeometry, elemVolVars, scv);
@@ -320,12 +320,12 @@ public:
      * \note has to be implemented by the model specific residual class
      *
      */
-    ResidualVector computeFlux(const Problem& problem,
-                               const Element& element,
-                               const FVElementGeometry& fvGeometry,
-                               const ElementVolumeVariables& elemVolVars,
-                               const SubControlVolumeFace& scvf,
-                               const ElementFluxVariablesCache& elemFluxVarsCache) const
+    NumEqVector computeFlux(const Problem& problem,
+                            const Element& element,
+                            const FVElementGeometry& fvGeometry,
+                            const ElementVolumeVariables& elemVolVars,
+                            const SubControlVolumeFace& scvf,
+                            const ElementFluxVariablesCache& elemFluxVarsCache) const
     {
         DUNE_THROW(Dune::NotImplemented, "This model does not implement a flux method!");
     }
@@ -372,8 +372,8 @@ public:
         // doing the time discretization...
 
         //! Compute storage with the model specific storage residual
-        ResidualVector prevStorage = asImp().computeStorage(problem, scv, prevVolVars);
-        ResidualVector storage = asImp().computeStorage(problem, scv, curVolVars);
+        NumEqVector prevStorage = asImp().computeStorage(problem, scv, prevVolVars);
+        NumEqVector storage = asImp().computeStorage(problem, scv, curVolVars);
 
         prevStorage *= prevVolVars.extrusionFactor();
         storage *= curVolVars.extrusionFactor();
@@ -407,7 +407,7 @@ public:
     {
         //! Compute source with the model specific storage residual
         const auto& curVolVars = curElemVolVars[scv];
-        ResidualVector source = asImp().computeSource(problem, element, fvGeometry, curElemVolVars, scv);
+        NumEqVector source = asImp().computeSource(problem, element, fvGeometry, curElemVolVars, scv);
         source *= scv.volume()*curVolVars.extrusionFactor();
 
         //! subtract source from local rate (sign convention in user interface)
@@ -452,12 +452,12 @@ public:
      * \param elemFluxVarsCache The flux variable caches for the element stencil
      * \param scvf The sub control volume face the flux term is integrated over
      */
-    ResidualVector evalFlux(const Problem& problem,
-                            const Element& element,
-                            const FVElementGeometry& fvGeometry,
-                            const ElementVolumeVariables& elemVolVars,
-                            const ElementFluxVariablesCache& elemFluxVarsCache,
-                            const SubControlVolumeFace& scvf) const
+    NumEqVector evalFlux(const Problem& problem,
+                         const Element& element,
+                         const FVElementGeometry& fvGeometry,
+                         const ElementVolumeVariables& elemVolVars,
+                         const ElementFluxVariablesCache& elemFluxVarsCache,
+                         const SubControlVolumeFace& scvf) const
     {
         return asImp().evalFlux(problem, element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
     }
diff --git a/dumux/common/fvproblem.hh b/dumux/common/fvproblem.hh
index be49c0f35192b8f6fa97d446470bd82c5359bac9..08e75df775acece6f2c14998e08333fff160e1a1 100644
--- a/dumux/common/fvproblem.hh
+++ b/dumux/common/fvproblem.hh
@@ -53,7 +53,7 @@ class FVProblem
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
@@ -248,10 +248,10 @@ public:
      * Negative values mean influx.
      * E.g. for the mass balance that would the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
      */
-    ResidualVector neumann(const Element& element,
-                           const FVElementGeometry& fvGeometry,
-                           const ElementVolumeVariables& elemVolVars,
-                           const SubControlVolumeFace& scvf) const
+    NumEqVector neumann(const Element& element,
+                        const FVElementGeometry& fvGeometry,
+                        const ElementVolumeVariables& elemVolVars,
+                        const SubControlVolumeFace& scvf) const
     {
         // forward it to the interface with only the global position
         return asImp_().neumannAtPos(scvf.ipGlobal());
@@ -266,11 +266,11 @@ public:
      * Negative values mean influx.
      * E.g. for the mass balance that would be the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
      */
-    ResidualVector neumannAtPos(const GlobalPosition &globalPos) const
+    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
     {
         //! As a default, i.e. if the user's problem does not overload any neumann method
         //! return no-flow Neumann boundary conditions at all Neumann boundaries
-        return ResidualVector(0.0);
+        return NumEqVector(0.0);
     }
 
     /*!
@@ -291,10 +291,10 @@ public:
      * that the conserved quantity is created, negative ones mean that it vanishes.
      * E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$.
      */
-    ResidualVector source(const Element &element,
-                          const FVElementGeometry& fvGeometry,
-                          const ElementVolumeVariables& elemVolVars,
-                          const SubControlVolume &scv) const
+    NumEqVector source(const Element &element,
+                       const FVElementGeometry& fvGeometry,
+                       const ElementVolumeVariables& elemVolVars,
+                       const SubControlVolume &scv) const
     {
         // forward to solution independent, fully-implicit specific interface
         return asImp_().sourceAtPos(scv.center());
@@ -313,11 +313,11 @@ public:
      * that the conserved quantity is created, negative ones mean that it vanishes.
      * E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$.
      */
-    ResidualVector sourceAtPos(const GlobalPosition &globalPos) const
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
     {
         //! As a default, i.e. if the user's problem does not overload any source method
         //! return 0.0 (no source terms)
-        return ResidualVector(0.0);
+        return NumEqVector(0.0);
     }
 
     /*!
@@ -399,12 +399,12 @@ public:
      *        Caution: Only overload this method in the implementation if you know
      *                 what you are doing.
      */
-    ResidualVector scvPointSources(const Element &element,
-                                   const FVElementGeometry& fvGeometry,
-                                   const ElementVolumeVariables& elemVolVars,
-                                   const SubControlVolume &scv) const
+    NumEqVector scvPointSources(const Element &element,
+                                const FVElementGeometry& fvGeometry,
+                                const ElementVolumeVariables& elemVolVars,
+                                const SubControlVolume &scv) const
     {
-        ResidualVector source(0);
+        NumEqVector source(0);
         auto scvIdx = scv.indexInElement();
         auto key = std::make_pair(fvGridGeometry_->elementMapper().index(element), scvIdx);
         if (pointSourceMap_.count(key))
diff --git a/dumux/common/pointsource.hh b/dumux/common/pointsource.hh
index 4faaf639c3e16b08ee5f85e0d4003da5b47e48df..7d5204ff5053f348a925d0513a0e190d2a361849 100644
--- a/dumux/common/pointsource.hh
+++ b/dumux/common/pointsource.hh
@@ -52,7 +52,7 @@ class PointSource
 {
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
@@ -64,7 +64,7 @@ class PointSource
 
 public:
     //! Constructor for constant point sources
-    PointSource(GlobalPosition pos, PrimaryVariables values)
+    PointSource(GlobalPosition pos, NumEqVector values)
       : values_(values), pos_(pos), embeddings_(1) {}
 
     //! Constructor for sol dependent point sources, when there is no
@@ -101,7 +101,7 @@ public:
     }
 
     //! Convenience = operator overload modifying only the values
-    PointSource& operator= (const PrimaryVariables& values)
+    PointSource& operator= (const NumEqVector& values)
     {
         values_ = values;
         return *this;
@@ -115,7 +115,7 @@ public:
     }
 
     //! return the source values
-    PrimaryVariables values() const
+    NumEqVector values() const
     { return values_; }
 
     //! return the source position
@@ -155,7 +155,7 @@ public:
     }
 
 protected:
-    PrimaryVariables values_; //!< value of the point source for each equation
+    NumEqVector values_; //!< value of the point source for each equation
 private:
     GlobalPosition pos_; //!< position of the point source
     std::size_t embeddings_; //!< how many SCVs the point source is associated with
@@ -171,27 +171,27 @@ class IdPointSource : public PointSource<TypeTag>
     using ParentType = PointSource<TypeTag>;
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
 
     static const int dimworld = GridView::dimensionworld;
     using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
 
 public:
     //! Constructor for constant point sources
-    IdPointSource(GlobalPosition pos, PrimaryVariables values, IdType id)
+    IdPointSource(GlobalPosition pos, NumEqVector values, IdType id)
       :  ParentType(pos, values), id_(id) {}
 
     //! Constructor for sol dependent point sources, when there is no
     // value known at the time of initialization
     IdPointSource(GlobalPosition pos, IdType id)
-      : ParentType(pos, PrimaryVariables(0.0)), id_(id) {}
+      : ParentType(pos, NumEqVector(0.0)), id_(id) {}
 
     //! return the sources identifier
     IdType id() const
     { return id_; }
 
     //! Convenience = operator overload modifying only the values
-    IdPointSource& operator= (const PrimaryVariables& values)
+    IdPointSource& operator= (const NumEqVector& values)
     {
         ParentType::operator=(values);
         return *this;
@@ -218,7 +218,7 @@ class SolDependentPointSource : public PointSource<TypeTag>
     using ParentType = PointSource<TypeTag>;
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
@@ -227,19 +227,19 @@ class SolDependentPointSource : public PointSource<TypeTag>
 
     static const int dimworld = GridView::dimensionworld;
     using GlobalPosition = typename Dune::FieldVector<Scalar, dimworld>;
-    // returns the PointSource values as PrimaryVariables
-    using ValueFunction = typename std::function<PrimaryVariables(const Problem &problem,
-                                                                  const Element &element,
-                                                                  const FVElementGeometry &fvGeometry,
-                                                                  const ElementVolumeVariables &elemVolVars,
-                                                                  const SubControlVolume &scv)>;
+    // returns the PointSource values as NumEqVector
+    using ValueFunction = typename std::function<NumEqVector(const Problem &problem,
+                                                             const Element &element,
+                                                             const FVElementGeometry &fvGeometry,
+                                                             const ElementVolumeVariables &elemVolVars,
+                                                             const SubControlVolume &scv)>;
 
 public:
     //! Constructor for sol dependent point sources, when there is no
     // value known at the time of initialization
     SolDependentPointSource(GlobalPosition pos,
                             ValueFunction valueFunction)
-      : ParentType(pos, PrimaryVariables(0.0)), valueFunction_(valueFunction) {}
+      : ParentType(pos, NumEqVector(0.0)), valueFunction_(valueFunction) {}
 
     //! an update function called before adding the value
     // to the local residual in the problem in scvPointSources
@@ -252,7 +252,7 @@ public:
     { this->values_ = valueFunction_(problem, element, fvGeometry, elemVolVars, scv); }
 
     //! Convenience = operator overload modifying only the values
-    SolDependentPointSource& operator= (const PrimaryVariables& values)
+    SolDependentPointSource& operator= (const NumEqVector& values)
     {
         ParentType::operator=(values);
         return *this;
diff --git a/dumux/porousmediumflow/2p1c/localresidual.hh b/dumux/porousmediumflow/2p1c/localresidual.hh
index 7a892ae01773498dd53cc9185dc7991162379dfa..6bbebc4dc31204b78aec644a504630eebe4d3ca0 100644
--- a/dumux/porousmediumflow/2p1c/localresidual.hh
+++ b/dumux/porousmediumflow/2p1c/localresidual.hh
@@ -39,7 +39,7 @@ class TwoPOneCLocalResidual : public ImmiscibleLocalResidual<TypeTag>
 {
     using ParentType = ImmiscibleLocalResidual<TypeTag>;
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
@@ -62,11 +62,11 @@ public:
 
 
     //! Evaluate the storage term within a given scv
-    ResidualVector computeStorage(const Problem& problem,
-                                  const SubControlVolume& scv,
-                                  const VolumeVariables& volVars) const
+    NumEqVector computeStorage(const Problem& problem,
+                               const SubControlVolume& scv,
+                               const VolumeVariables& volVars) const
     {
-        ResidualVector storage(0.0);
+        NumEqVector storage(0.0);
         // Compute storage term of all components within all phases
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
         {
@@ -84,17 +84,17 @@ public:
     }
 
     //! Evaluate the fluxes over a face of a sub control volume
-    ResidualVector computeFlux(const Problem& problem,
-                               const Element& element,
-                               const FVElementGeometry& fvGeometry,
-                               const ElementVolumeVariables& elemVolVars,
-                               const SubControlVolumeFace& scvf,
-                               const ElementFluxVariablesCache& elemFluxVarsCache) const
+    NumEqVector computeFlux(const Problem& problem,
+                            const Element& element,
+                            const FVElementGeometry& fvGeometry,
+                            const ElementVolumeVariables& elemVolVars,
+                            const SubControlVolumeFace& scvf,
+                            const ElementFluxVariablesCache& elemFluxVarsCache) const
     {
         FluxVariables fluxVars;
         fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
 
-        ResidualVector flux;
+        NumEqVector flux;
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
         {
             // The physical quantities for which we perform upwinding
diff --git a/dumux/porousmediumflow/3p3c/localresidual.hh b/dumux/porousmediumflow/3p3c/localresidual.hh
index b3c76110a0d0dd5a205ae241fff49f1c8db0d44d..df127346e46c423ca8ad7240478f3d053d4d7e15 100644
--- a/dumux/porousmediumflow/3p3c/localresidual.hh
+++ b/dumux/porousmediumflow/3p3c/localresidual.hh
@@ -46,7 +46,7 @@ class ThreePThreeCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
     using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
     using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
@@ -88,11 +88,11 @@ public:
      * \param scv The sub control volume
      * \param volVars The volume variables
      */
-    ResidualVector computeStorage(const Problem& problem,
-                                  const SubControlVolume& scv,
-                                  const VolumeVariables& volVars) const
+    NumEqVector computeStorage(const Problem& problem,
+                               const SubControlVolume& scv,
+                               const VolumeVariables& volVars) const
     {
-        ResidualVector storage(0.0);
+        NumEqVector storage(0.0);
 
         // compute storage term of all components within all phases
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
@@ -127,18 +127,18 @@ public:
      * \param scvf The sub control volume face
      * \param elemFluxVarsCache The element flux variables cache
      */
-    ResidualVector computeFlux(const Problem& problem,
-                               const Element& element,
-                               const FVElementGeometry& fvGeometry,
-                               const ElementVolumeVariables& elemVolVars,
-                               const SubControlVolumeFace& scvf,
-                               const ElementFluxVariablesCache& elemFluxVarsCache) const
+    NumEqVector computeFlux(const Problem& problem,
+                            const Element& element,
+                            const FVElementGeometry& fvGeometry,
+                            const ElementVolumeVariables& elemVolVars,
+                            const SubControlVolumeFace& scvf,
+                            const ElementFluxVariablesCache& elemFluxVarsCache) const
     {
         FluxVariables fluxVars;
         fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
 
         // get upwind weights into local scope
-        ResidualVector flux(0.0);
+        NumEqVector flux(0.0);
 
         // advective fluxes
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
diff --git a/dumux/porousmediumflow/3pwateroil/localresidual.hh b/dumux/porousmediumflow/3pwateroil/localresidual.hh
index 2a6e504db3d906b07b7fb09254cf43a5b2f2d2d3..d4f456af651419b56dc38bcc6306f8095c28330f 100644
--- a/dumux/porousmediumflow/3pwateroil/localresidual.hh
+++ b/dumux/porousmediumflow/3pwateroil/localresidual.hh
@@ -48,7 +48,7 @@ protected:
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
     using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
     using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
@@ -90,11 +90,11 @@ public:
      *  \param scvIdx The SCV (sub-control-volume) index
      *  \param usePrevSol Evaluate function with solution of current or previous time step
      */
-     ResidualVector computeStorage(const Problem& problem,
-                                  const SubControlVolume& scv,
-                                  const VolumeVariables& volVars) const
+     NumEqVector computeStorage(const Problem& problem,
+                                const SubControlVolume& scv,
+                                const VolumeVariables& volVars) const
     {
-        ResidualVector storage(0.0);
+        NumEqVector storage(0.0);
 
         const auto massOrMoleDensity = [](const auto& volVars, const int phaseIdx)
         { return useMoles ? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); };
@@ -133,18 +133,18 @@ public:
      * \param onBoundary A boolean variable to specify whether the flux variables
      *        are calculated for interior SCV faces or boundary faces, default=false
      */
-    ResidualVector computeFlux(const Problem& problem,
-                               const Element& element,
-                               const FVElementGeometry& fvGeometry,
-                               const ElementVolumeVariables& elemVolVars,
-                               const SubControlVolumeFace& scvf,
-                               const ElementFluxVariablesCache& elemFluxVarsCache) const
+    NumEqVector computeFlux(const Problem& problem,
+                            const Element& element,
+                            const FVElementGeometry& fvGeometry,
+                            const ElementVolumeVariables& elemVolVars,
+                            const SubControlVolumeFace& scvf,
+                            const ElementFluxVariablesCache& elemFluxVarsCache) const
     {
         FluxVariables fluxVars;
         fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
 
         // get upwind weights into local scope
-        ResidualVector flux(0.0);
+        NumEqVector flux(0.0);
 
         const auto massOrMoleDensity = [](const auto& volVars, const int phaseIdx)
         { return useMoles ? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); };
diff --git a/dumux/porousmediumflow/compositional/localresidual.hh b/dumux/porousmediumflow/compositional/localresidual.hh
index dcb303ec835a254a5c465c92ea9f97b85ce57dee..098e70f9dfff7af529cc39f2097f3050e877a244 100644
--- a/dumux/porousmediumflow/compositional/localresidual.hh
+++ b/dumux/porousmediumflow/compositional/localresidual.hh
@@ -45,7 +45,7 @@ class CompositionalLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidua
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
     using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
     using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
@@ -81,11 +81,11 @@ public:
      * \param scv The sub control volume
      * \param volVars The current or previous volVars
      */
-    ResidualVector computeStorage(const Problem& problem,
-                                  const SubControlVolume& scv,
-                                  const VolumeVariables& volVars) const
+    NumEqVector computeStorage(const Problem& problem,
+                               const SubControlVolume& scv,
+                               const VolumeVariables& volVars) const
     {
-        ResidualVector storage(0.0);
+        NumEqVector storage(0.0);
 
         const auto massOrMoleDensity = [](const auto& volVars, const int phaseIdx)
         { return useMoles ? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); };
@@ -133,17 +133,17 @@ public:
      * \param scvf The sub control volume face to compute the flux on
      * \param elemFluxVarsCache The cache related to flux compuation
      */
-    ResidualVector computeFlux(const Problem& problem,
-                               const Element& element,
-                               const FVElementGeometry& fvGeometry,
-                               const ElementVolumeVariables& elemVolVars,
-                               const SubControlVolumeFace& scvf,
-                               const ElementFluxVariablesCache& elemFluxVarsCache) const
+    NumEqVector computeFlux(const Problem& problem,
+                            const Element& element,
+                            const FVElementGeometry& fvGeometry,
+                            const ElementVolumeVariables& elemVolVars,
+                            const SubControlVolumeFace& scvf,
+                            const ElementFluxVariablesCache& elemFluxVarsCache) const
     {
         FluxVariables fluxVars;
         fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
         // get upwind weights into local scope
-        ResidualVector flux(0.0);
+        NumEqVector flux(0.0);
 
         const auto massOrMoleDensity = [](const auto& volVars, const int phaseIdx)
         { return useMoles ? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); };
diff --git a/dumux/porousmediumflow/immiscible/localresidual.hh b/dumux/porousmediumflow/immiscible/localresidual.hh
index 200079cc4f7b60b7e8f37b03f8c5c6ab10a495de..0b60c613ff730b6c03a8107241fe893415256d4a 100644
--- a/dumux/porousmediumflow/immiscible/localresidual.hh
+++ b/dumux/porousmediumflow/immiscible/localresidual.hh
@@ -40,7 +40,7 @@ class ImmiscibleLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
     using ParentType = typename GET_PROP_TYPE(TypeTag, BaseLocalResidual);
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
@@ -71,12 +71,12 @@ public:
      * \note The volVars can be different to allow computing
      *       the implicit euler time derivative here
      */
-    ResidualVector computeStorage(const Problem& problem,
-                                  const SubControlVolume& scv,
-                                  const VolumeVariables& volVars) const
+    NumEqVector computeStorage(const Problem& problem,
+                               const SubControlVolume& scv,
+                               const VolumeVariables& volVars) const
     {
         // partial time derivative of the phase mass
-        ResidualVector storage;
+        NumEqVector storage;
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
         {
             auto eqIdx = conti0EqIdx + phaseIdx;
@@ -105,17 +105,17 @@ public:
      * \param scvf The sub control volume face to compute the flux on
      * \param elemFluxVarsCache The cache related to flux compuation
      */
-    ResidualVector computeFlux(const Problem& problem,
-                               const Element& element,
-                               const FVElementGeometry& fvGeometry,
-                               const ElementVolumeVariables& elemVolVars,
-                               const SubControlVolumeFace& scvf,
-                               const ElementFluxVariablesCache& elemFluxVarsCache) const
+    NumEqVector computeFlux(const Problem& problem,
+                            const Element& element,
+                            const FVElementGeometry& fvGeometry,
+                            const ElementVolumeVariables& elemVolVars,
+                            const SubControlVolumeFace& scvf,
+                            const ElementFluxVariablesCache& elemFluxVarsCache) const
     {
         FluxVariables fluxVars;
         fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
 
-        ResidualVector flux;
+        NumEqVector flux;
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
         {
             // the physical quantities for which we perform upwinding
diff --git a/dumux/porousmediumflow/mineralization/localresidual.hh b/dumux/porousmediumflow/mineralization/localresidual.hh
index 0eddfe8366e930de1e403a9c1d7350f23000b3a0..f99d8af5a731a12d274607a2d47b9ae8d4f70927 100644
--- a/dumux/porousmediumflow/mineralization/localresidual.hh
+++ b/dumux/porousmediumflow/mineralization/localresidual.hh
@@ -38,7 +38,7 @@ template<class TypeTag>
 class MineralizationLocalResidual: public CompositionalLocalResidual<TypeTag>
 {
     using ParentType = CompositionalLocalResidual<TypeTag>;
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
@@ -70,9 +70,9 @@ public:
      * \param volVars The volume variables (primary/secondary variables) in the scv
      * \return Amount per volume of the conserved quantities
      */
-    ResidualVector computeStorage(const Problem& problem,
-                                  const SubControlVolume& scv,
-                                  const VolumeVariables& volVars) const
+    NumEqVector computeStorage(const Problem& problem,
+                               const SubControlVolume& scv,
+                               const VolumeVariables& volVars) const
     {
         auto storage = ParentType::computeStorage(problem, scv, volVars);
 
diff --git a/dumux/porousmediumflow/nonequilibrium/localresidual.hh b/dumux/porousmediumflow/nonequilibrium/localresidual.hh
index 20d59a876ff4a7c5ddd6c64be439461b7cb502c4..87031244378919f58e186c3d0315a250c2f2dd94 100644
--- a/dumux/porousmediumflow/nonequilibrium/localresidual.hh
+++ b/dumux/porousmediumflow/nonequilibrium/localresidual.hh
@@ -47,7 +47,7 @@ class NonEquilibriumLocalResidualImplementation<TypeTag, true, false>: public GE
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using ParentType = typename GET_PROP_TYPE(TypeTag, EquilibriumLocalResidual);
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
@@ -74,17 +74,17 @@ public:
      *
      */
 
-    ResidualVector computeFlux(const Problem& problem,
-                               const Element& element,
-                               const FVElementGeometry& fvGeometry,
-                               const ElementVolumeVariables& elemVolVars,
-                               const SubControlVolumeFace& scvf,
-                               const ElementFluxVariablesCache& elemFluxVarsCache) const
+    NumEqVector computeFlux(const Problem& problem,
+                            const Element& element,
+                            const FVElementGeometry& fvGeometry,
+                            const ElementVolumeVariables& elemVolVars,
+                            const SubControlVolumeFace& scvf,
+                            const ElementFluxVariablesCache& elemFluxVarsCache) const
     {
         FluxVariables fluxVars;
         fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
         // get upwind weights into local scope
-        ResidualVector flux(0.0);
+        NumEqVector flux(0.0);
 
         const auto moleDensity = [](const auto& volVars, const int phaseIdx)
         { return volVars.molarDensity(phaseIdx); };
@@ -122,13 +122,13 @@ public:
 
     }
 
-    ResidualVector computeSource(const Problem& problem,
-                                 const Element& element,
-                                 const FVElementGeometry& fvGeometry,
-                                 const ElementVolumeVariables& elemVolVars,
-                                 const SubControlVolume &scv) const
+    NumEqVector computeSource(const Problem& problem,
+                              const Element& element,
+                              const FVElementGeometry& fvGeometry,
+                              const ElementVolumeVariables& elemVolVars,
+                              const SubControlVolume &scv) const
     {
-        ResidualVector source(0.0);
+        NumEqVector source(0.0);
 
         // add contributions from volume flux sources
         source += problem.source(element, fvGeometry, elemVolVars, scv);
@@ -161,7 +161,7 @@ class NonEquilibriumLocalResidualImplementation<TypeTag, true, true>: public GET
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
     using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
     using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
@@ -194,11 +194,11 @@ public:
      *    \param storage The mass of the component within the sub-control volume
      *    \param volVars The volume variables
      */
-    ResidualVector computeStorage(const Problem& problem,
-                                  const SubControlVolume& scv,
-                                  const VolumeVariables& volVars) const
+    NumEqVector computeStorage(const Problem& problem,
+                               const SubControlVolume& scv,
+                               const VolumeVariables& volVars) const
     {
-       ResidualVector storage(0.0);
+       NumEqVector storage(0.0);
 
      // compute storage term of all components within all phases
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
@@ -226,17 +226,17 @@ public:
      *        \param fluxVars The flux Variables
      *        \param elemVolVars The volume variables of the current element
      */
-    ResidualVector computeFlux(const Problem& problem,
-                               const Element& element,
-                               const FVElementGeometry& fvGeometry,
-                               const ElementVolumeVariables& elemVolVars,
-                               const SubControlVolumeFace& scvf,
-                               const ElementFluxVariablesCache& elemFluxVarsCache) const
+    NumEqVector computeFlux(const Problem& problem,
+                            const Element& element,
+                            const FVElementGeometry& fvGeometry,
+                            const ElementVolumeVariables& elemVolVars,
+                            const SubControlVolumeFace& scvf,
+                            const ElementFluxVariablesCache& elemFluxVarsCache) const
     {
         FluxVariables fluxVars;
         fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
         // get upwind weights into local scope
-        ResidualVector flux(0.0);
+        NumEqVector flux(0.0);
 
         // advective fluxes
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
@@ -268,13 +268,13 @@ public:
     }
 
 
-    ResidualVector computeSource(const Problem& problem,
-                                 const Element& element,
-                                 const FVElementGeometry& fvGeometry,
-                                 const ElementVolumeVariables& elemVolVars,
-                                 const SubControlVolume &scv) const
+    NumEqVector computeSource(const Problem& problem,
+                              const Element& element,
+                              const FVElementGeometry& fvGeometry,
+                              const ElementVolumeVariables& elemVolVars,
+                              const SubControlVolume &scv) const
     {
-         ResidualVector source(0.0);
+         NumEqVector source(0.0);
         // In the case of a kinetic consideration, mass transfer
         // between phases is realized via source terms there is a
         // balance equation for each component in each phase
diff --git a/dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh b/dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh
index 0b92e94363df8b9601cb063c52abe108d780956c..ee8c2976177b7636d755b7460f5c22dc40b908da 100644
--- a/dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh
+++ b/dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh
@@ -44,7 +44,7 @@ template<class TypeTag>
 class EnergyLocalResidualNonEquilibrium<TypeTag, 1/*numEnergyEqFluid*/>
 {
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
@@ -70,7 +70,7 @@ public:
 
 
     //! The energy storage in the fluid phase with index phaseIdx
-    static void fluidPhaseStorage(ResidualVector& storage,
+    static void fluidPhaseStorage(NumEqVector& storage,
                                   const SubControlVolume& scv,
                                   const VolumeVariables& volVars,
                                   int phaseIdx)
@@ -85,7 +85,7 @@ public:
 
 
     //! The energy storage in the solid matrix
-    static void solidPhaseStorage(ResidualVector& storage,
+    static void solidPhaseStorage(NumEqVector& storage,
                                   const SubControlVolume& scv,
                                   const VolumeVariables& volVars)
     {
@@ -100,14 +100,14 @@ public:
     }
 
      // this is to make nonequilibrium work with compositional local residual, compositional calls that for non-isothermal models
-    static void heatConvectionFlux(ResidualVector& flux,
+    static void heatConvectionFlux(NumEqVector& flux,
                                    FluxVariables& fluxVars,
                                    int phaseIdx)
     {}
 
 
    //! The advective phase energy fluxes
-    static void heatConvectionFlux(ResidualVector& flux,
+    static void heatConvectionFlux(NumEqVector& flux,
                                    FluxVariables& fluxVars,
                                    const ElementVolumeVariables& elemVolVars,
                                    const SubControlVolumeFace& scvf,
@@ -142,7 +142,7 @@ public:
     }
 
     //! The diffusive energy fluxes
-    static void heatConductionFlux(ResidualVector& flux,
+    static void heatConductionFlux(NumEqVector& flux,
                                    FluxVariables& fluxVars)
     {
         //in case we have one energy equation for more than one fluid phase we use an effective law in the nonequilibrium fourierslaw
@@ -157,7 +157,7 @@ public:
      *
      * \param scv The sub-control volume over which we integrate the source term
      */
-    static void computeSourceEnergy(ResidualVector& source,
+    static void computeSourceEnergy(NumEqVector& source,
                                     const Element& element,
                                     const FVElementGeometry& fvGeometry,
                                     const ElementVolumeVariables& elemVolVars,
@@ -313,7 +313,7 @@ class EnergyLocalResidualNonEquilibrium<TypeTag, 2 /*numEnergyEqFluid*/>
 : public EnergyLocalResidualNonEquilibrium<TypeTag, 1 /*numEnergyEqFluid*/>
 {
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
@@ -342,7 +342,7 @@ class EnergyLocalResidualNonEquilibrium<TypeTag, 2 /*numEnergyEqFluid*/>
 public:
 
     //! The energy storage in the fluid phase with index phaseIdx
-    static void fluidPhaseStorage(ResidualVector& storage,
+    static void fluidPhaseStorage(NumEqVector& storage,
                                   const SubControlVolume& scv,
                                   const VolumeVariables& volVars,
                                   int phaseIdx)
@@ -355,7 +355,7 @@ public:
 
 
    //! The advective phase energy fluxes
-    static void heatConvectionFlux(ResidualVector& flux,
+    static void heatConvectionFlux(NumEqVector& flux,
                                    FluxVariables& fluxVars,
                                    const ElementVolumeVariables& elemVolVars,
                                    const SubControlVolumeFace& scvf,
@@ -390,7 +390,7 @@ public:
     }
 
     //! The diffusive energy fluxes
-    static void heatConductionFlux(ResidualVector& flux,
+    static void heatConductionFlux(NumEqVector& flux,
                                    FluxVariables& fluxVars)
     {
 
@@ -406,7 +406,7 @@ public:
      *
      * \param scv The sub-control volume over which we integrate the source term
      */
-    static void computeSourceEnergy(ResidualVector& source,
+    static void computeSourceEnergy(NumEqVector& source,
                                     const Element& element,
                                     const FVElementGeometry& fvGeometry,
                                     const ElementVolumeVariables& elemVolVars,
diff --git a/dumux/porousmediumflow/nonisothermal/localresidual.hh b/dumux/porousmediumflow/nonisothermal/localresidual.hh
index cba5152c46859780eb8c7ae2ec4c0382f6d8faa1..115d3569aefc9a5446658af81fa85ab25124ed4f 100644
--- a/dumux/porousmediumflow/nonisothermal/localresidual.hh
+++ b/dumux/porousmediumflow/nonisothermal/localresidual.hh
@@ -44,7 +44,7 @@ template<class TypeTag>
 class EnergyLocalResidualImplementation<TypeTag, false>
 {
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
@@ -59,7 +59,7 @@ public:
      * \param volVars The volume variables
      * \param phaseIdx The phase index
      */
-    static void fluidPhaseStorage(ResidualVector& storage,
+    static void fluidPhaseStorage(NumEqVector& storage,
                                   const SubControlVolume& scv,
                                   const VolumeVariables& volVars,
                                   int phaseIdx)
@@ -72,7 +72,7 @@ public:
      * \param scv The sub-control volume
      * \param volVars The volume variables
      */
-    static void solidPhaseStorage(ResidualVector& storage,
+    static void solidPhaseStorage(NumEqVector& storage,
                                   const SubControlVolume& scv,
                                   const VolumeVariables& volVars)
     {}
@@ -84,7 +84,7 @@ public:
      * \param fluxVars The flux variables.
      * \param phaseIdx The phase index
      */
-    static void heatConvectionFlux(ResidualVector& flux,
+    static void heatConvectionFlux(NumEqVector& flux,
                                    FluxVariables& fluxVars,
                                    int phaseIdx)
     {}
@@ -95,7 +95,7 @@ public:
      * \param flux TODO docme!
      * \param fluxVars The flux variables.
      */
-    static void heatConductionFlux(ResidualVector& flux,
+    static void heatConductionFlux(NumEqVector& flux,
                                    FluxVariables& fluxVars)
     {}
 };
@@ -108,7 +108,7 @@ template<class TypeTag>
 class EnergyLocalResidualImplementation<TypeTag, true>
 {
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
@@ -127,7 +127,7 @@ public:
      * \param volVars The volume variables
      * \param phaseIdx The phase index
      */
-    static void fluidPhaseStorage(ResidualVector& storage,
+    static void fluidPhaseStorage(NumEqVector& storage,
                                   const SubControlVolume& scv,
                                   const VolumeVariables& volVars,
                                   int phaseIdx)
@@ -145,7 +145,7 @@ public:
      * \param scv The sub-control volume
      * \param volVars The volume variables
      */
-    static void solidPhaseStorage(ResidualVector& storage,
+    static void solidPhaseStorage(NumEqVector& storage,
                                   const SubControlVolume& scv,
                                   const VolumeVariables& volVars)
     {
@@ -162,7 +162,7 @@ public:
      * \param fluxVars The flux variables.
      * \param phaseIdx The phase index
      */
-    static void heatConvectionFlux(ResidualVector& flux,
+    static void heatConvectionFlux(NumEqVector& flux,
                                    FluxVariables& fluxVars,
                                    int phaseIdx)
     {
@@ -178,7 +178,7 @@ public:
      * \param flux TODO docme!
      * \param fluxVars The flux variables.
      */
-    static void heatConductionFlux(ResidualVector& flux,
+    static void heatConductionFlux(NumEqVector& flux,
                                    FluxVariables& fluxVars)
     {
         flux[energyEqIdx] += fluxVars.heatConductionFlux();
diff --git a/dumux/porousmediumflow/richards/localresidual.hh b/dumux/porousmediumflow/richards/localresidual.hh
index 85f2d4af63703525e7cbad5705b67cc6e1ce989f..a74420e0d2e8f7b5f87de24c74383220bfad9e95 100644
--- a/dumux/porousmediumflow/richards/localresidual.hh
+++ b/dumux/porousmediumflow/richards/localresidual.hh
@@ -42,7 +42,7 @@ class RichardsLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
     using ParentType = typename GET_PROP_TYPE(TypeTag, BaseLocalResidual);
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
@@ -81,12 +81,12 @@ public:
      * \note The volVars can be different to allow computing
      *       the implicit euler time derivative here
      */
-    ResidualVector computeStorage(const Problem& problem,
-                                  const SubControlVolume& scv,
-                                  const VolumeVariables& volVars) const
+    NumEqVector computeStorage(const Problem& problem,
+                               const SubControlVolume& scv,
+                               const VolumeVariables& volVars) const
     {
         // partial time derivative of the phase mass
-        ResidualVector storage(0.0);
+        NumEqVector storage(0.0);
         storage[conti0EqIdx] = volVars.porosity()
                                * volVars.density(wPhaseIdx)
                                * volVars.saturation(wPhaseIdx);
@@ -118,17 +118,17 @@ public:
      * \param scvf The sub control volume face to compute the flux on
      * \param elemFluxVarsCache The cache related to flux compuation
      */
-    ResidualVector computeFlux(const Problem& problem,
-                               const Element& element,
-                               const FVElementGeometry& fvGeometry,
-                               const ElementVolumeVariables& elemVolVars,
-                               const SubControlVolumeFace& scvf,
-                               const ElementFluxVariablesCache& elemFluxVarsCache) const
+    NumEqVector computeFlux(const Problem& problem,
+                            const Element& element,
+                            const FVElementGeometry& fvGeometry,
+                            const ElementVolumeVariables& elemVolVars,
+                            const SubControlVolumeFace& scvf,
+                            const ElementFluxVariablesCache& elemFluxVarsCache) const
     {
         FluxVariables fluxVars;
         fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
 
-        ResidualVector flux(0.0);
+        NumEqVector flux(0.0);
         // the physical quantities for which we perform upwinding
         auto upwindTerm = [](const auto& volVars)
                           { return volVars.density(wPhaseIdx)*volVars.mobility(wPhaseIdx); };
diff --git a/dumux/porousmediumflow/tracer/localresidual.hh b/dumux/porousmediumflow/tracer/localresidual.hh
index ea1fa8d1bcb7d5cf5fed7be8fe9dd56c1ad15e80..ea845f1c59bad54e849c1e346a2b6166a50fa015 100644
--- a/dumux/porousmediumflow/tracer/localresidual.hh
+++ b/dumux/porousmediumflow/tracer/localresidual.hh
@@ -47,7 +47,7 @@ class TracerLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
-    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
     using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
@@ -74,11 +74,11 @@ public:
      * \param scv The sub control volume
      * \param volVars The primary and secondary varaibles on the scv
      */
-    PrimaryVariables computeStorage(const Problem& problem,
-                                    const SubControlVolume& scv,
-                                    const VolumeVariables& volVars) const
+    NumEqVector computeStorage(const Problem& problem,
+                               const SubControlVolume& scv,
+                               const VolumeVariables& volVars) const
     {
-        PrimaryVariables storage(0.0);
+        NumEqVector storage(0.0);
 
         // formulation with mole balances
         if (useMoles)
@@ -111,17 +111,17 @@ public:
      * \param scvf The sub control volume face
      * \param elemFluxVarsCache The cache related to flux compuation
      */
-    PrimaryVariables computeFlux(const Problem& problem,
-                                 const Element& element,
-                                 const FVElementGeometry& fvGeometry,
-                                 const ElementVolumeVariables& elemVolVars,
-                                 const SubControlVolumeFace& scvf,
-                                 const ElementFluxVariablesCache& elemFluxVarsCache) const
+    NumEqVector computeFlux(const Problem& problem,
+                            const Element& element,
+                            const FVElementGeometry& fvGeometry,
+                            const ElementVolumeVariables& elemVolVars,
+                            const SubControlVolumeFace& scvf,
+                            const ElementFluxVariablesCache& elemFluxVarsCache) const
     {
         FluxVariables fluxVars;
         fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
 
-        PrimaryVariables flux(0.0);
+        NumEqVector flux(0.0);
         const auto diffusiveFluxes = fluxVars.molecularDiffusionFlux(phaseIdx);
 
         // formulation with mole balances