diff --git a/test/freeflow/navierstokes/angelitestproblem.hh b/test/freeflow/navierstokes/angelitestproblem.hh
index a4d1305c666d8e8bfae5bd2db0ec668cf293ed2b..a6561feb94ae6a47a1f24260be959a9963320b6c 100644
--- a/test/freeflow/navierstokes/angelitestproblem.hh
+++ b/test/freeflow/navierstokes/angelitestproblem.hh
@@ -98,7 +98,7 @@ class AngeliTestProblem : public NavierStokesProblem<TypeTag>
     using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using SourceValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
     using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
 
@@ -161,9 +161,9 @@ public:
      *
      * \param globalPos The global position
      */
-    SourceValues sourceAtPos(const GlobalPosition &globalPos) const
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
     {
-        return SourceValues(0.0);
+        return NumEqVector(0.0);
     }
 
     // \}
diff --git a/test/freeflow/navierstokes/channeltestproblem.hh b/test/freeflow/navierstokes/channeltestproblem.hh
index bc893ef8cb21c0c219aceaad9c3a9dbfdbeb6d86..8dfe135e62116d050512e9705d134daed5dcd849 100644
--- a/test/freeflow/navierstokes/channeltestproblem.hh
+++ b/test/freeflow/navierstokes/channeltestproblem.hh
@@ -110,7 +110,7 @@ class ChannelTestProblem : public NavierStokesProblem<TypeTag>
     using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using SourceValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
 
     using TimeLoopPtr = std::shared_ptr<CheckPointTimeLoop<Scalar>>;
 
@@ -145,9 +145,9 @@ public:
      *
      * \param globalPos The global position
      */
-    SourceValues sourceAtPos(const GlobalPosition &globalPos) const
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
     {
-        return SourceValues(0.0);
+        return NumEqVector(0.0);
     }
     // \}
    /*!
diff --git a/test/freeflow/navierstokes/closedsystemtestproblem.hh b/test/freeflow/navierstokes/closedsystemtestproblem.hh
index 35c455f70e434cdd952e6fc4b6017f85a170ea4f..3ecf2fcbe8c48bebb6b5c8cb8389cbffb59b7503 100644
--- a/test/freeflow/navierstokes/closedsystemtestproblem.hh
+++ b/test/freeflow/navierstokes/closedsystemtestproblem.hh
@@ -91,7 +91,7 @@ class ClosedSystemTestProblem : public NavierStokesProblem<TypeTag>
     using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using SourceValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
 
 public:
@@ -130,9 +130,9 @@ public:
      * \param values Stores the source values, acts as return value
      * \param globalPos The global position
      */
-    SourceValues sourceAtPos(const GlobalPosition &globalPos) const
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
     {
-        return SourceValues(0.0);
+        return NumEqVector(0.0);
     }
     // \}
    /*!
diff --git a/test/freeflow/navierstokes/doneatestproblem.hh b/test/freeflow/navierstokes/doneatestproblem.hh
index 1613a194de8984a89b084f571bf8b73affd86b9f..5ba0d15881c54bb8f6ca57cd5bdad2c1f83073f3 100644
--- a/test/freeflow/navierstokes/doneatestproblem.hh
+++ b/test/freeflow/navierstokes/doneatestproblem.hh
@@ -103,7 +103,7 @@ class DoneaTestProblem : public NavierStokesProblem<TypeTag>
     using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using SourceValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
 
     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
 
@@ -172,9 +172,9 @@ public:
      *
      * \param globalPos The global position
      */
-    SourceValues sourceAtPos(const GlobalPosition &globalPos) const
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
     {
-        SourceValues source(0.0);
+        NumEqVector source(0.0);
         Scalar x = globalPos[0];
         Scalar y = globalPos[1];
 
diff --git a/test/freeflow/navierstokes/kovasznaytestproblem.hh b/test/freeflow/navierstokes/kovasznaytestproblem.hh
index 27f3925dd81b3cfddfbbfe7ac2ba1f432a147aa4..ce99d609fdf58e17c81c65c1270443f2f1fca9e0 100644
--- a/test/freeflow/navierstokes/kovasznaytestproblem.hh
+++ b/test/freeflow/navierstokes/kovasznaytestproblem.hh
@@ -96,7 +96,7 @@ class KovasznayTestProblem : public NavierStokesProblem<TypeTag>
     using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using SourceValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
 
     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
 
@@ -164,9 +164,9 @@ public:
      *
      * \param globalPos The global position
      */
-    SourceValues sourceAtPos(const GlobalPosition &globalPos) const
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
     {
-        return SourceValues(0.0);
+        return NumEqVector(0.0);
     }
 
     // \}
diff --git a/test/freeflow/navierstokesnc/channeltestproblem.hh b/test/freeflow/navierstokesnc/channeltestproblem.hh
index f999ee6f8e6c9697908b8b621fbbf1fae84d8119..5cebba5386f4cd52ef2a27306866bc32b725f008 100644
--- a/test/freeflow/navierstokesnc/channeltestproblem.hh
+++ b/test/freeflow/navierstokesnc/channeltestproblem.hh
@@ -122,7 +122,7 @@ class ChannelNCTestProblem : public NavierStokesProblem<TypeTag>
     using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using SourceValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
 
     using TimeLoopPtr = std::shared_ptr<CheckPointTimeLoop<Scalar>>;
     using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
@@ -161,9 +161,9 @@ public:
      *
      * \param globalPos The global position
      */
-    SourceValues sourceAtPos(const GlobalPosition &globalPos) const
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
     {
-        return SourceValues(0.0);
+        return NumEqVector(0.0);
     }
     // \}
    /*!
diff --git a/test/freeflow/navierstokesnc/densityflowproblem.hh b/test/freeflow/navierstokesnc/densityflowproblem.hh
index 3af036a487b05175aa95a7ecadec0cad01773d1e..bd6fc5e41094dd863ef5f68efd186d87b81368c7 100644
--- a/test/freeflow/navierstokesnc/densityflowproblem.hh
+++ b/test/freeflow/navierstokesnc/densityflowproblem.hh
@@ -111,7 +111,7 @@ class DensityDrivenFlowProblem : public NavierStokesProblem<TypeTag>
     using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using SourceValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
 
     using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
@@ -152,9 +152,9 @@ public:
      *
      * \param globalPos The global position
      */
-    SourceValues sourceAtPos(const GlobalPosition &globalPos) const
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
     {
-        return SourceValues(0.0);
+        return NumEqVector(0.0);
     }
     // \}
    /*!
diff --git a/test/freeflow/navierstokesnc/msfreeflowtestproblem.hh b/test/freeflow/navierstokesnc/msfreeflowtestproblem.hh
index 9f55cb50b241729627299a2da5f6b14da201b518..9f008b1fb086a80063b3b62c752d32ced8424f92 100644
--- a/test/freeflow/navierstokesnc/msfreeflowtestproblem.hh
+++ b/test/freeflow/navierstokesnc/msfreeflowtestproblem.hh
@@ -207,7 +207,7 @@ class MaxwellStefanNCTestProblem : public NavierStokesProblem<TypeTag>
     using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using SourceValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
 
     using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
@@ -352,9 +352,9 @@ public:
      *
      * \param globalPos The global position
      */
-    SourceValues sourceAtPos(const GlobalPosition &globalPos) const
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
     {
-        return SourceValues(0.0);
+        return NumEqVector(0.0);
     }
 
     // \}
diff --git a/test/porousmediumflow/1p/implicit/1pniconvectionproblem.hh b/test/porousmediumflow/1p/implicit/1pniconvectionproblem.hh
index 8356ddee48b8ad5dd52eb33f0c963719e8c2b093..898a40d7962cff658ed5ad731f3b599d75da8f32 100644
--- a/test/porousmediumflow/1p/implicit/1pniconvectionproblem.hh
+++ b/test/porousmediumflow/1p/implicit/1pniconvectionproblem.hh
@@ -125,7 +125,7 @@ class OnePNIConvectionProblem : public PorousMediumFlowProblem<TypeTag>
         energyEqIdx = Indices::energyEqIdx
     };
 
-    using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using Element = typename GridView::template Codim<0>::Entity;
     using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
@@ -269,12 +269,12 @@ public:
      * The \a values store the mass flux of each phase normal to the boundary.
      * Negative values indicate an inflow.
      */
-    NeumannFluxes neumann(const Element& element,
+    NumEqVector neumann(const Element& element,
                              const FVElementGeometry& fvGeometry,
                              const ElementVolumeVariables& elemVolvars,
                              const SubControlVolumeFace& scvf) const
     {
-        NeumannFluxes values(0.0);
+        NumEqVector values(0.0);
         const auto globalPos = scvf.ipGlobal();
         const auto& volVars = elemVolvars[scvf.insideScvIdx()];
 
diff --git a/test/porousmediumflow/1p/implicit/tubesproblem.hh b/test/porousmediumflow/1p/implicit/tubesproblem.hh
index 85c835e47188489e0a5909f825edb9334c733c85..919ddc44f9b51140494f17962c4f0539a3112315 100644
--- a/test/porousmediumflow/1p/implicit/tubesproblem.hh
+++ b/test/porousmediumflow/1p/implicit/tubesproblem.hh
@@ -93,6 +93,7 @@ class TubesTestProblem : public PorousMediumFlowProblem<TypeTag>
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using Element = typename GridView::template Codim<0>::Entity;
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
@@ -210,12 +211,12 @@ 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$.
      */
-    PrimaryVariables source(const Element &element,
+    NumEqVector source(const Element &element,
                             const FVElementGeometry& fvGeometry,
                             const ElementVolumeVariables& elemVolVars,
                             const SubControlVolume &scv) const
     {
-        PrimaryVariables source(0.0);
+        NumEqVector source(0.0);
 
         const auto& globalPos = scv.center();
         const auto& volVars = elemVolVars[scv];
diff --git a/test/porousmediumflow/1pnc/implicit/1p2cniconductionproblem.hh b/test/porousmediumflow/1pnc/implicit/1p2cniconductionproblem.hh
index 6a3cc069c552a27ac72ae65e619504c1388860fe..550a5bcf9e7c98ae58156a3e55b333b68a8f3f7a 100644
--- a/test/porousmediumflow/1pnc/implicit/1p2cniconductionproblem.hh
+++ b/test/porousmediumflow/1pnc/implicit/1p2cniconductionproblem.hh
@@ -107,7 +107,7 @@ class OnePTwoCNIConductionProblem : public PorousMediumFlowProblem<TypeTag>
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
     using Element = typename GridView::template Codim<0>::Entity;
     using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);
@@ -259,8 +259,8 @@ public:
      * \brief Evaluate the boundary conditions for a neumann
      *        boundary segment.
      */
-    ResidualVector neumannAtPos(const GlobalPosition &globalPos) const
-    { return ResidualVector(0.0); }
+    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
+    { return NumEqVector(0.0); }
 
     // \}
 
@@ -280,8 +280,8 @@ public:
      *
      * The units must be according to either using mole or mass fractions. (mole/(m^3*s) or kg/(m^3*s))
      */
-    PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
-    { return PrimaryVariables(0.0); }
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+    { return NumEqVector(0.0); }
 
     /*!
      * \brief Evaluate the initial value for a control volume.
diff --git a/test/porousmediumflow/1pnc/implicit/1p2cniconvectionproblem.hh b/test/porousmediumflow/1pnc/implicit/1p2cniconvectionproblem.hh
index 475fcec9f57877a53ca836dcdf8505827b732ea7..4b9e76c1f4391238ac1e8b21b1e68a70ffc906e7 100644
--- a/test/porousmediumflow/1pnc/implicit/1p2cniconvectionproblem.hh
+++ b/test/porousmediumflow/1pnc/implicit/1p2cniconvectionproblem.hh
@@ -108,7 +108,7 @@ class OnePTwoCNIConvectionProblem : public PorousMediumFlowProblem<TypeTag>
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
     using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
@@ -277,12 +277,12 @@ public:
      * in normal direction of each phase. 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,
+    NumEqVector neumann(const Element& element,
                            const FVElementGeometry& fvGeometry,
                            const ElementVolumeVariables& elemVolVars,
                            const SubControlVolumeFace& scvf) const
     {
-        ResidualVector flux(0.0);
+        NumEqVector flux(0.0);
         const auto& globalPos = scvf.ipGlobal();
         const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
 
@@ -315,8 +315,8 @@ public:
      *
      * The units must be according to either using mole or mass fractions. (mole/(m^3*s) or kg/(m^3*s))
      */
-    PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
-    { return PrimaryVariables(0.0); }
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+    { return NumEqVector(0.0); }
 
     /*!
      * \brief Evaluate the initial value for a control volume.
diff --git a/test/porousmediumflow/1pnc/implicit/1p2ctestproblem.hh b/test/porousmediumflow/1pnc/implicit/1p2ctestproblem.hh
index 555f3e24193187e6b44f9c6e476efad1a9ec74df..d7349db4d17002996eb221608dd965ad80f357c7 100644
--- a/test/porousmediumflow/1pnc/implicit/1p2ctestproblem.hh
+++ b/test/porousmediumflow/1pnc/implicit/1p2ctestproblem.hh
@@ -109,7 +109,7 @@ class OnePTwoCTestProblem : public PorousMediumFlowProblem<TypeTag>
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
     using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
@@ -225,7 +225,7 @@ public:
      * in normal direction of each phase. 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,
+    NumEqVector neumann(const Element& element,
                            const FVElementGeometry& fvGeometry,
                            const ElementVolumeVariables& elemVolVars,
                            const SubControlVolumeFace& scvf) const
@@ -233,7 +233,7 @@ public:
         // set a fixed pressure on the right side of the domain
         const Scalar dirichletPressure = 1e5;
 
-        ResidualVector flux(0.0);
+        NumEqVector flux(0.0);
         const auto& ipGlobal = scvf.ipGlobal();
         const auto& volVars = elemVolVars[scvf.insideScvIdx()];
 
@@ -317,8 +317,8 @@ public:
      *
      * The units must be according to either using mole or mass fractions. (mole/(m^3*s) or kg/(m^3*s))
      */
-    PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
-    { return PrimaryVariables(0.0); }
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+    { return NumEqVector(0.0); }
 
     /*!
      * \brief Evaluate the initial value for a control volume.
diff --git a/test/porousmediumflow/1pncmin/implicit/thermochemproblem.hh b/test/porousmediumflow/1pncmin/implicit/thermochemproblem.hh
index c72e474f658c7a7ef308d019f17df9445fa72b9c..7bfbe7838065406537f0e23d4a90160ca4058c59 100644
--- a/test/porousmediumflow/1pncmin/implicit/thermochemproblem.hh
+++ b/test/porousmediumflow/1pncmin/implicit/thermochemproblem.hh
@@ -97,7 +97,7 @@ class ThermoChemProblem : public PorousMediumFlowProblem<TypeTag>
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
     using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
     using ReactionRate =ThermoChemReaction<TypeTag>;
@@ -218,12 +218,12 @@ public:
      * Negative values indicate an inflow.
      */
 
-    ResidualVector neumann(const Element& element,
+    NumEqVector neumann(const Element& element,
                            const FVElementGeometry& fvGeometry,
                            const ElementVolumeVariables& elemVolVars,
                            const SubControlVolumeFace& scvf) const
     {
-        ResidualVector flux(0.0);
+        NumEqVector flux(0.0);
         return flux;
     }
 
@@ -280,13 +280,13 @@ 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$.
      */
-    PrimaryVariables source(const Element &element,
+    NumEqVector source(const Element &element,
                             const FVElementGeometry& fvGeometry,
                             const ElementVolumeVariables& elemVolVars,
                             const SubControlVolume &scv) const
     {
 
-        PrimaryVariables source(0.0);
+        NumEqVector source(0.0);
         const auto& volVars = elemVolVars[scv];
 
         Scalar qMass = rrate_.thermoChemReaction(volVars);
diff --git a/test/porousmediumflow/2p/implicit/fracture/problem.hh b/test/porousmediumflow/2p/implicit/fracture/problem.hh
index 7dab629e6f8fd3132b4cb4e0941839c2ff4e9d94..859ea96a764daa62e22c3a0cc1cf0d94cd17073f 100644
--- a/test/porousmediumflow/2p/implicit/fracture/problem.hh
+++ b/test/porousmediumflow/2p/implicit/fracture/problem.hh
@@ -90,6 +90,7 @@ class FractureProblem : public PorousMediumFlowProblem<TypeTag>
     using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
 
     enum {
 
@@ -147,9 +148,9 @@ public:
      *               \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$
      * \param globalPos The global position
      */
-    PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
     {
-        return PrimaryVariables(0.0);
+        return NumEqVector(0.0);
     }
 
     // \}
@@ -208,9 +209,9 @@ public:
      * For this method, the \a values parameter stores the mass flux
      * in normal direction of each phase. Negative values mean influx.
      */
-    PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
+    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
     {
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
         if (onInlet_(globalPos)) {
             values[contiNEqIdx] = -0.04; // kg / (m * s)
         }
diff --git a/test/porousmediumflow/2p/implicit/incompressible/problem.hh b/test/porousmediumflow/2p/implicit/incompressible/problem.hh
index 5abe543119a0676e029b1c7b2aedbc14839b5707..a3de10baaae3202b03b43d75b23d5cf16455d9c6 100644
--- a/test/porousmediumflow/2p/implicit/incompressible/problem.hh
+++ b/test/porousmediumflow/2p/implicit/incompressible/problem.hh
@@ -90,7 +90,7 @@ class TwoPTestProblem : public PorousMediumFlowProblem<TypeTag>
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
     using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimensionworld>;
-    using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
     enum {
         pwIdx = Indices::pwIdx,
@@ -161,9 +161,9 @@ public:
      * For this method, the \a values parameter stores the mass flux
      * in normal direction of each phase. Negative values mean influx.
      */
-    NeumannFluxes neumannAtPos(const GlobalPosition &globalPos) const
+    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
     {
-        NeumannFluxes values(0.0);
+        NumEqVector values(0.0);
         if (onInlet_(globalPos))
             values[contiNEqIdx] = -0.04; // kg / (m * s)
         return values;
diff --git a/test/porousmediumflow/2p/implicit/nonisothermal/problem.hh b/test/porousmediumflow/2p/implicit/nonisothermal/problem.hh
index 0fc40e9725cbbc82af3824971fc8a0f266155e58..f35f3dd1f740700a73c25e1ebbbff4b4e41b40f7 100644
--- a/test/porousmediumflow/2p/implicit/nonisothermal/problem.hh
+++ b/test/porousmediumflow/2p/implicit/nonisothermal/problem.hh
@@ -121,8 +121,7 @@ class InjectionProblem2PNI : public PorousMediumFlowProblem<TypeTag>
     };
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
-    using Sources = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
 
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
@@ -176,9 +175,9 @@ public:
      *               \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$
      * \param globalPos The global position
      */
-    Sources sourceAtPos(const GlobalPosition &globalPos) const
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
     {
-        Sources values(0.0);
+        NumEqVector values(0.0);
         values = 0;
         return values;
     }
@@ -242,12 +241,12 @@ public:
      * The \a values store the mass flux of each phase normal to the boundary.
      * Negative values indicate an inflow.
      */
-    NeumannFluxes neumann(const Element &element,
+    NumEqVector neumann(const Element &element,
                           const FVElementGeometry& fvGeometry,
                           const ElementVolumeVariables& elemVolVars,
                           const SubControlVolumeFace& scvf) const
     {
-        NeumannFluxes values(0.0);
+        NumEqVector values(0.0);
         const auto globalPos = scvf.ipGlobal();
 
         if (globalPos[1] < 13.75 + eps_ && globalPos[1] > 6.875 - eps_)
diff --git a/test/porousmediumflow/2p1c/implicit/steaminjectionproblem.hh b/test/porousmediumflow/2p1c/implicit/steaminjectionproblem.hh
index b8b515e9d5633330a6f3e7b8d11ca3c55428ad56..124a0f56e9bef640411ed90f9a11e84a082e0f17 100644
--- a/test/porousmediumflow/2p1c/implicit/steaminjectionproblem.hh
+++ b/test/porousmediumflow/2p1c/implicit/steaminjectionproblem.hh
@@ -86,7 +86,7 @@ class InjectionProblem : public PorousMediumFlowProblem<TypeTag>
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using Sources = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
@@ -94,7 +94,6 @@ class InjectionProblem : public PorousMediumFlowProblem<TypeTag>
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
     using Element = typename GridView::template Codim<0>::Entity;
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
-    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
 
     // copy some indices for convenience
     enum {
@@ -130,12 +129,12 @@ public:
 
 
     //! \copydoc Dumux::FVProblem::source()
-    Sources source(const Element &element,
+    NumEqVector source(const Element &element,
                    const FVElementGeometry& fvGeometry,
                    const ElementVolumeVariables& elemVolVars,
                    const SubControlVolume &scv) const
     {
-          return Sources(0.0);
+          return NumEqVector(0.0);
     }
 
     /*!
@@ -191,12 +190,12 @@ public:
      * in normal direction of each phase. 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,
+    NumEqVector neumann(const Element& element,
                            const FVElementGeometry& fvGeometry,
                            const ElementVolumeVariables& elemVolVars,
                            const SubControlVolumeFace& scvf) const
     {
-        ResidualVector values(0.0);
+        NumEqVector values(0.0);
 
         const auto& ipGlobal = scvf.ipGlobal();
 
diff --git a/test/porousmediumflow/2p2c/implicit/injectionproblem.hh b/test/porousmediumflow/2p2c/implicit/injectionproblem.hh
index a1420e0e560b9cb253ff9b41b99a176192c41b93..e39f4765ae794bb77ca45261f3b7ffe2b5206411 100644
--- a/test/porousmediumflow/2p2c/implicit/injectionproblem.hh
+++ b/test/porousmediumflow/2p2c/implicit/injectionproblem.hh
@@ -121,7 +121,7 @@ class InjectionProblem : public PorousMediumFlowProblem<TypeTag>
     };
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
     using Element = typename GridView::template Codim<0>::Entity;
@@ -241,12 +241,12 @@ public:
      * The \a values store the mass flux of each phase normal to the boundary.
      * Negative values indicate an inflow.
      */
-    NeumannFluxes neumann(const Element& element,
+    NumEqVector neumann(const Element& element,
                           const FVElementGeometry& fvGeometry,
                           const ElementVolumeVariables& elemVolVars,
                           const SubControlVolumeFace& scvf) const
     {
-        NeumannFluxes values(0.0);
+        NumEqVector values(0.0);
 
         const auto& globalPos = scvf.ipGlobal();
 
diff --git a/test/porousmediumflow/2p2c/implicit/waterairproblem.hh b/test/porousmediumflow/2p2c/implicit/waterairproblem.hh
index 60797abd7ca220635ba1b8c47ef0ecf6e4c5d851..1fb629d56705aa3a89590ee9ba413d42f0b02cca 100644
--- a/test/porousmediumflow/2p2c/implicit/waterairproblem.hh
+++ b/test/porousmediumflow/2p2c/implicit/waterairproblem.hh
@@ -128,7 +128,7 @@ class WaterAirProblem : public PorousMediumFlowProblem<TypeTag>
     };
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
     using Element = typename GridView::template Codim<0>::Entity;
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
@@ -232,12 +232,12 @@ public:
      *
      * The units must be according to either using mole or mass fractions. (mole/(m^2*s) or kg/(m^2*s))
      */
-    NeumannFluxes neumann(const Element& element,
+    NumEqVector neumann(const Element& element,
                           const FVElementGeometry& fvGeometry,
                           const ElementVolumeVariables& elemVolVars,
                           const SubControlVolumeFace& scvf) const
     {
-        NeumannFluxes values(0.0);
+        NumEqVector values(0.0);
 
         const auto& globalPos = scvf.ipGlobal();
 
diff --git a/test/porousmediumflow/2pnc/implicit/2pncdiffusionproblem.hh b/test/porousmediumflow/2pnc/implicit/2pncdiffusionproblem.hh
index 3539d15ebf72947aadfa3a301e6499bfd8e265b8..33de974b751ee203669a9d572e915d4736aa4fee 100644
--- a/test/porousmediumflow/2pnc/implicit/2pncdiffusionproblem.hh
+++ b/test/porousmediumflow/2pnc/implicit/2pncdiffusionproblem.hh
@@ -107,7 +107,7 @@ class TwoPNCDiffusionProblem : public PorousMediumFlowProblem<TypeTag>
     };
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
     using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
@@ -212,9 +212,9 @@ public:
      *
      * The units must be according to either using mole or mass fractions. (mole/(m^2*s) or kg/(m^2*s))
      */
-    NeumannFluxes neumannAtPos(const GlobalPosition& globalPos) const
+    NumEqVector neumannAtPos(const GlobalPosition& globalPos) const
     {
-        NeumannFluxes values(0.0);
+        NumEqVector values(0.0);
         return values;
     }
 
diff --git a/test/porousmediumflow/2pnc/implicit/fuelcellproblem.hh b/test/porousmediumflow/2pnc/implicit/fuelcellproblem.hh
index b0312a144d3d56dba9db334f3220892f4e7c2196..71b18bad20a1489afa27d18ba2998327000cc436 100644
--- a/test/porousmediumflow/2pnc/implicit/fuelcellproblem.hh
+++ b/test/porousmediumflow/2pnc/implicit/fuelcellproblem.hh
@@ -84,7 +84,7 @@ class FuelCellProblem : public PorousMediumFlowProblem<TypeTag>
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using Sources = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
@@ -167,12 +167,12 @@ public:
     { return temperature_; }
 
     //! \copydoc Dumux::FVProblem::source()
-    Sources source(const Element &element,
+    NumEqVector source(const Element &element,
                    const FVElementGeometry& fvGeometry,
                    const ElementVolumeVariables& elemVolVars,
                    const SubControlVolume &scv) const
     {
-        Sources values(0.0);
+        NumEqVector values(0.0);
         const auto& globalPos = scv.dofPosition();
 
         //reaction sources from electro chemistry
diff --git a/test/porousmediumflow/2pncmin/implicit/dissolutionproblem.hh b/test/porousmediumflow/2pncmin/implicit/dissolutionproblem.hh
index 9a720af2b9cdd4742f03932bb2bf1c1991f9aa13..0e83673507d309a94eb798d870a7b21a844e7d0d 100644
--- a/test/porousmediumflow/2pncmin/implicit/dissolutionproblem.hh
+++ b/test/porousmediumflow/2pncmin/implicit/dissolutionproblem.hh
@@ -125,7 +125,7 @@ class DissolutionProblem : public PorousMediumFlowProblem<TypeTag>
     };
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using Sources = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using Element = typename GridView::template Codim<0>::Entity;
@@ -314,12 +314,12 @@ 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$.
      */
-    Sources source(const Element &element,
+    NumEqVector source(const Element &element,
                    const FVElementGeometry& fvGeometry,
                    const ElementVolumeVariables& elemVolVars,
                    const SubControlVolume &scv) const
     {
-        Sources source(0.0);
+        NumEqVector source(0.0);
 
         const auto& volVars = elemVolVars[scv];
 
diff --git a/test/porousmediumflow/3p/implicit/3pniconductionproblem.hh b/test/porousmediumflow/3p/implicit/3pniconductionproblem.hh
index d93cd2d90b85fed61ac5b38cf0137e4dbafaf712..648ee3c0449c2ee622fea6f21c05a96c56d16d14 100644
--- a/test/porousmediumflow/3p/implicit/3pniconductionproblem.hh
+++ b/test/porousmediumflow/3p/implicit/3pniconductionproblem.hh
@@ -116,7 +116,7 @@ class ThreePNIConductionProblem : public PorousMediumFlowProblem<TypeTag>
     using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
     using IapwsH2O = H2O<Scalar>;
-    using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
 
     // copy some indices for convenience
     using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
@@ -262,9 +262,9 @@ public:
      * For this method, the \a values parameter stores the mass flux
      * in normal direction of each phase. Negative values mean influx.
      */
-    NeumannFluxes neumannAtPos(const GlobalPosition &globalPos) const
+    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
     {
-        return NeumannFluxes(0.0);
+        return NumEqVector(0.0);
     }
 
     // \}
@@ -286,9 +286,9 @@ public:
      *
      * The units must be according to either using mole or mass fractions. (mole/(m^3*s) or kg/(m^3*s))
      */
-    PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
     {
-        return PrimaryVariables(0.0);
+        return NumEqVector(0.0);
     }
 
     /*!
diff --git a/test/porousmediumflow/3p/implicit/3pniconvectionproblem.hh b/test/porousmediumflow/3p/implicit/3pniconvectionproblem.hh
index 67a0383bcb406a9a1cfe1f5a9a1588414ff7b283..05d9674f396a08eaba17221040fd1d6bb5cf6d04 100644
--- a/test/porousmediumflow/3p/implicit/3pniconvectionproblem.hh
+++ b/test/porousmediumflow/3p/implicit/3pniconvectionproblem.hh
@@ -108,6 +108,7 @@ class ThreePNIConvectionProblem : public PorousMediumFlowProblem<TypeTag>
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
     using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
@@ -266,12 +267,12 @@ public:
      * \param scvf The subcontrolvolume face
      *  Negative values mean influx.
      */
-    PrimaryVariables neumann(const Element &element,
+    NumEqVector neumann(const Element &element,
                              const FVElementGeometry& fvGeometry,
                              const ElementVolumeVariables& elemVolVars,
                              const SubControlVolumeFace& scvf) const
     {
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
         const auto globalPos = scvf.ipGlobal();
         const auto& volVars = elemVolVars[scvf.insideScvIdx()];
 
diff --git a/test/porousmediumflow/3p/implicit/infiltration3pproblem.hh b/test/porousmediumflow/3p/implicit/infiltration3pproblem.hh
index c7f7e9a8ea0dc118b78d5b4ffdca384c13c8406f..86e8ad3d56a6eb7e1ef084902e011ef97a2d230b 100644
--- a/test/porousmediumflow/3p/implicit/infiltration3pproblem.hh
+++ b/test/porousmediumflow/3p/implicit/infiltration3pproblem.hh
@@ -113,7 +113,7 @@ class InfiltrationThreePProblem : public PorousMediumFlowProblem<TypeTag>
     };
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
@@ -222,9 +222,9 @@ public:
      * For this method, the \a values parameter stores the mass flux
      * in normal direction of each phase. Negative values mean influx.
      */
-    NeumannFluxes neumannAtPos(const GlobalPosition &globalPos) const
+    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
     {
-        NeumannFluxes values(0.0);
+        NumEqVector values(0.0);
 
         // negative values for injection
         if (time_ < 2592000.0 - eps_)
diff --git a/test/porousmediumflow/3p3c/implicit/columnxylolproblem.hh b/test/porousmediumflow/3p3c/implicit/columnxylolproblem.hh
index cf209daa1d78aebba636feeb8d957333b374070e..6b266e6163625873228b682c24804481a4a26205 100644
--- a/test/porousmediumflow/3p3c/implicit/columnxylolproblem.hh
+++ b/test/porousmediumflow/3p3c/implicit/columnxylolproblem.hh
@@ -114,7 +114,7 @@ class ColumnProblem : public PorousMediumFlowProblem<TypeTag>
     };
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
     using Element = typename GridView::template Codim<0>::Entity;
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
@@ -202,12 +202,12 @@ public:
      * For this method, the \a values parameter stores the mass flux
      * in normal direction of each phase. Negative values mean influx.
      */
-    NeumannFluxes neumann(const Element& element,
+    NumEqVector neumann(const Element& element,
                           const FVElementGeometry& fvGeometry,
                           const ElementVolumeVariables& elemVolVars,
                           const SubControlVolumeFace& scvf) const
     {
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
         const auto& globalPos = scvf.ipGlobal();
 
         // negative values for injection
diff --git a/test/porousmediumflow/3p3c/implicit/infiltration3p3cproblem.hh b/test/porousmediumflow/3p3c/implicit/infiltration3p3cproblem.hh
index 2ac4f5212d9f39fda116763571063af6716c5eb1..e3394126f162200095402946871afe322bcfc882 100644
--- a/test/porousmediumflow/3p3c/implicit/infiltration3p3cproblem.hh
+++ b/test/porousmediumflow/3p3c/implicit/infiltration3p3cproblem.hh
@@ -120,7 +120,7 @@ class InfiltrationThreePThreeCProblem : public PorousMediumFlowProblem<TypeTag>
     };
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
@@ -223,9 +223,9 @@ public:
      * For this method, the \a values parameter stores the mass flux
      * in normal direction of each phase. Negative values mean influx.
      */
-    NeumannFluxes neumannAtPos(const GlobalPosition &globalPos) const
+    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
     {
-        NeumannFluxes values(0.0);
+        NumEqVector values(0.0);
 
         // negative values for injection
         if ((globalPos[0] < 80.0 + eps_) && (globalPos[0] > 55.0 - eps_) && (globalPos[1] > 10.0 - eps_))
diff --git a/test/porousmediumflow/3p3c/implicit/kuevetteproblem.hh b/test/porousmediumflow/3p3c/implicit/kuevetteproblem.hh
index 63d53d14e72883c33a1b8c4a1c2292703a5ff6c1..eca71a09304dc72818303c5e6a8c96e93645bda9 100644
--- a/test/porousmediumflow/3p3c/implicit/kuevetteproblem.hh
+++ b/test/porousmediumflow/3p3c/implicit/kuevetteproblem.hh
@@ -130,7 +130,7 @@ class KuevetteProblem : public PorousMediumFlowProblem<TypeTag>
     };
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
     using Element = typename GridView::template Codim<0>::Entity;
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
@@ -219,12 +219,12 @@ public:
      * For this method, the \a values parameter stores the mass flux
      * in normal direction of each phase. Negative values mean influx.
      */
-    NeumannFluxes neumann(const Element& element,
+    NumEqVector neumann(const Element& element,
                           const FVElementGeometry& fvGeometry,
                           const ElementVolumeVariables& elemVolVars,
                           const SubControlVolumeFace& scvf) const
     {
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
         const auto& globalPos = scvf.ipGlobal();
 
         // negative values for injection
diff --git a/test/porousmediumflow/3pwateroil/implicit/3pwateroilsagdproblem.hh b/test/porousmediumflow/3pwateroil/implicit/3pwateroilsagdproblem.hh
index 4c5ebb1f1eed245b88fe53e3c591da1cfb20e742..8fd14728f0697459536a1d9cd1e4ca50bb3fdf8b 100644
--- a/test/porousmediumflow/3pwateroil/implicit/3pwateroilsagdproblem.hh
+++ b/test/porousmediumflow/3pwateroil/implicit/3pwateroilsagdproblem.hh
@@ -101,7 +101,7 @@ class SagdProblem : public PorousMediumFlowProblem<TypeTag>
     };
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
     using Element = typename GridView::template Codim<0>::Entity;
@@ -216,12 +216,12 @@ public:
      * For this method, the \a values parameter stores the mass flux
      * in normal direction of each phase. Negative values mean influx.
      */
-    NeumannFluxes neumann(const Element& element,
+    NumEqVector neumann(const Element& element,
                           const FVElementGeometry& fvGeometry,
                           const ElementVolumeVariables& elemVolVars,
                           const SubControlVolumeFace& scvf) const
     {
-        NeumannFluxes values(0.0);
+        NumEqVector values(0.0);
 
          const auto& globalPos = scvf.ipGlobal();
         // negative values for injection at injection well
diff --git a/test/porousmediumflow/co2/implicit/heterogeneousproblem.hh b/test/porousmediumflow/co2/implicit/heterogeneousproblem.hh
index e5bd91d321a65647e15f7df6a6652a3f11835aba..a68bd010e91310e4fbee7f13ef63147c4abb48c0 100644
--- a/test/porousmediumflow/co2/implicit/heterogeneousproblem.hh
+++ b/test/porousmediumflow/co2/implicit/heterogeneousproblem.hh
@@ -153,7 +153,7 @@ class HeterogeneousProblem : public PorousMediumFlowProblem<TypeTag>
 #endif
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
     using Element = typename GridView::template Codim<0>::Entity;
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
@@ -363,14 +363,14 @@ public:
      * in normal direction of each phase. Negative values mean influx.
      * E.g. for the mass balance that would the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
      */
-    NeumannFluxes neumann(const Element& element,
+    NumEqVector neumann(const Element& element,
                           const FVElementGeometry& fvGeometry,
                           const ElementVolumeVariables& elemVolvars,
                           const SubControlVolumeFace& scvf) const
     {
         const auto boundaryId = scvf.boundaryFlag();
 
-        NeumannFluxes fluxes(0.0);
+        NumEqVector fluxes(0.0);
          // kg/(m^2*s) or mole/(m^2*s) depending on useMoles
         if (boundaryId == injectionBottom_)
         {
diff --git a/test/porousmediumflow/mpnc/implicit/combustionproblem1c.hh b/test/porousmediumflow/mpnc/implicit/combustionproblem1c.hh
index b67d06fe2339e04db3654f62459fc6dde7f8ba61..59a80669e08f9b6fd5758fc91326e1a986d9ca6e 100644
--- a/test/porousmediumflow/mpnc/implicit/combustionproblem1c.hh
+++ b/test/porousmediumflow/mpnc/implicit/combustionproblem1c.hh
@@ -117,6 +117,7 @@ class CombustionProblemOneComponent: public PorousMediumFlowProblem<TypeTag>
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
@@ -216,12 +217,12 @@ public:
      * Positive values mean that mass is created, negative ones mean that it vanishes.
      */
     //! \copydoc Dumux::ImplicitProblem::source()
-    PrimaryVariables source(const Element &element,
+    NumEqVector source(const Element &element,
                             const FVElementGeometry& fvGeometry,
                             const ElementVolumeVariables& elemVolVars,
                             const SubControlVolume &scv) const
     {
-        PrimaryVariables priVars(0.0);
+        NumEqVector values(0.0);
 
         const auto& globalPos = scv.dofPosition();
 
@@ -233,10 +234,10 @@ public:
             if (onRightBoundaryPorousMedium_(globalPos))
             {
                 // Testing the location of a vertex, but function is called for each associated scv. Compensate for that
-                priVars[energyEqSolidIdx] = heatFluxFromRight_ / volume / numScv;
+                values[energyEqSolidIdx] = heatFluxFromRight_ / volume / numScv;
             }
          }
-        return priVars;
+        return values;
     }
 
     /*!
@@ -287,12 +288,12 @@ public:
      * For this method, the \a values parameter stores the mass flux
      * in normal direction of each phase. Negative values mean influx.
      */
-    PrimaryVariables neumann(const Element &element,
+    NumEqVector neumann(const Element &element,
                              const FVElementGeometry& fvGeometry,
                              const ElementVolumeVariables& elemVolVars,
                              const SubControlVolumeFace& scvf) const
     {
-        PrimaryVariables priVars(0.0);
+        NumEqVector values(0.0);
 
         const auto& globalPos = fvGeometry.scv(scvf.insideScvIdx()).dofPosition();
         const auto& scvIdx = scvf.insideScvIdx();
@@ -327,11 +328,11 @@ public:
 
         if (onLeftBoundary_(globalPos))
         {
-            priVars[conti00EqIdx + wCompIdx] = - molarFlux * fluidState.moleFraction(wPhaseIdx, wCompIdx);;
-            priVars[conti00EqIdx + nCompIdx] = - molarFlux * fluidState.moleFraction(wPhaseIdx, nCompIdx);;
-            priVars[energyEq0Idx] = - massFluxInjectedPhase * fluidState.enthalpy(wPhaseIdx);
+            values[conti00EqIdx + wCompIdx] = - molarFlux * fluidState.moleFraction(wPhaseIdx, wCompIdx);;
+            values[conti00EqIdx + nCompIdx] = - molarFlux * fluidState.moleFraction(wPhaseIdx, nCompIdx);;
+            values[energyEq0Idx] = - massFluxInjectedPhase * fluidState.enthalpy(wPhaseIdx);
         }
-        return priVars;
+        return values;
     }
 
     /*!
diff --git a/test/porousmediumflow/mpnc/implicit/evaporationatmosphereproblem.hh b/test/porousmediumflow/mpnc/implicit/evaporationatmosphereproblem.hh
index 6205c76f098ebd45f3a0b3459a598955e07eaeae..25d711f61a709f05107405dd889fc5b26b98e5b1 100644
--- a/test/porousmediumflow/mpnc/implicit/evaporationatmosphereproblem.hh
+++ b/test/porousmediumflow/mpnc/implicit/evaporationatmosphereproblem.hh
@@ -105,6 +105,7 @@ class EvaporationAtmosphereProblem: public PorousMediumFlowProblem<TypeTag>
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
@@ -247,12 +248,12 @@ public:
      * The \a values store the mass flux of each phase normal to the boundary.
      * Negative values indicate an inflow.
      */
-    PrimaryVariables neumann(const Element &element,
+    NumEqVector neumann(const Element &element,
                              const FVElementGeometry& fvGeometry,
                              const ElementVolumeVariables& elemVolVars,
                              const SubControlVolumeFace& scvf) const
     {
-        PrimaryVariables priVars(0.0);
+        NumEqVector values(0.0);
         const auto& globalPos = fvGeometry.scv(scvf.insideScvIdx()).dofPosition();
         const Scalar massFluxInjectedPhase = massFluxInjectedPhase_ ;
 
@@ -293,15 +294,15 @@ public:
         // actually setting the fluxes
         if (onLeftBoundary_(globalPos) && this->spatialParams().inFF_(globalPos))
         {
-            priVars[conti00EqIdx + nPhaseIdx * numComponents + wCompIdx]
+            values[conti00EqIdx + nPhaseIdx * numComponents + wCompIdx]
              = -molarFlux * fluidState.moleFraction(nPhaseIdx, wCompIdx);
-            priVars[conti00EqIdx + nPhaseIdx * numComponents + nCompIdx]
+            values[conti00EqIdx + nPhaseIdx * numComponents + nCompIdx]
              = -molarFlux * fluidState.moleFraction(nPhaseIdx, nCompIdx);
             // energy equations are specified mass specifically
-            priVars[energyEq0Idx + nPhaseIdx] = - massFluxInjectedPhase
+            values[energyEq0Idx + nPhaseIdx] = - massFluxInjectedPhase
                                                     * fluidState.enthalpy(nPhaseIdx) ;
         }
-        return priVars;
+        return values;
     }
 
     /*!
@@ -333,13 +334,13 @@ public:
      *
      *      Positive values mean that mass is created, negative ones mean that it vanishes.
      */
-    PrimaryVariables source(const Element &element,
+    NumEqVector source(const Element &element,
                             const FVElementGeometry& fvGeometry,
                             const ElementVolumeVariables& elemVolVars,
                             const SubControlVolume &scv) const
     {
-        PrimaryVariables priVars(0.0);
-        return priVars;
+        NumEqVector values(0.0);
+        return values;
 
     }
 
diff --git a/test/porousmediumflow/mpnc/implicit/obstacleproblem.hh b/test/porousmediumflow/mpnc/implicit/obstacleproblem.hh
index 66a856afc9b4f9cb58bf4b90ba9c1f16ef0bb317..e7537ee9763e7bde0ebe4d1e961822ffb7e58d27 100644
--- a/test/porousmediumflow/mpnc/implicit/obstacleproblem.hh
+++ b/test/porousmediumflow/mpnc/implicit/obstacleproblem.hh
@@ -112,6 +112,7 @@ class ObstacleProblem
     using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
@@ -239,12 +240,12 @@ public:
      *
      * Negative values mean influx.
      */
-    PrimaryVariables neumann(const Element& element,
+    NumEqVector neumann(const Element& element,
                              const FVElementGeometry& fvGeometry,
                              const ElementVolumeVariables& elemVolVars,
                              const SubControlVolumeFace& scvf) const
     {
-        return PrimaryVariables(0.0);
+        return NumEqVector(0.0);
     }
 
     // \}
@@ -267,12 +268,12 @@ public:
      * Positive values mean that mass is created, negative ones mean that it vanishes.
      */
     //! \copydoc Dumux::ImplicitProblem::source()
-    PrimaryVariables source(const Element &element,
+    NumEqVector source(const Element &element,
                             const FVElementGeometry& fvGeometry,
                             const ElementVolumeVariables& elemVolVars,
                             const SubControlVolume &scv) const
     {
-       return PrimaryVariables(0.0);
+       return NumEqVector(0.0);
     }
 
     /*!
diff --git a/test/porousmediumflow/richards/implicit/richardsanalyticalproblem.hh b/test/porousmediumflow/richards/implicit/richardsanalyticalproblem.hh
index 1a83dbddfffe1b564e06c1ec4512a26e3ea6d22a..386e469ca006688cf68d950ce3099a284223f2ff 100644
--- a/test/porousmediumflow/richards/implicit/richardsanalyticalproblem.hh
+++ b/test/porousmediumflow/richards/implicit/richardsanalyticalproblem.hh
@@ -91,6 +91,7 @@ class RichardsAnalyticalProblem :  public PorousMediumFlowProblem<TypeTag>
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
@@ -172,9 +173,9 @@ public:
      * \param values Storage for all primary variables of the source term
      * \param globalPos The position for which the source term is set
      */
-    PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
     {
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
         const Scalar time = time_;
         const Scalar pwTop = 98942.8;
         const Scalar pwBottom = 95641.1;
@@ -256,9 +257,9 @@ public:
      * \param values The neumann values for the conservation equations
      * \param globalPos The position for which the Neumann value is set
      */
-    PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
+    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
     {
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
         return values;
     }
 
diff --git a/test/porousmediumflow/richards/implicit/richardslensproblem.hh b/test/porousmediumflow/richards/implicit/richardslensproblem.hh
index 1acf47f5f63e99f8b188bc128839577ecc9db899..7fd5a143019770103302e4d2df63c62b222f78a9 100644
--- a/test/porousmediumflow/richards/implicit/richardslensproblem.hh
+++ b/test/porousmediumflow/richards/implicit/richardslensproblem.hh
@@ -97,6 +97,7 @@ class RichardsLensProblem : public PorousMediumFlowProblem<TypeTag>
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
     using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
@@ -208,9 +209,9 @@ public:
      *
      * \param globalPos The position for which the Neumann value is set
      */
-    PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
+    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
     {
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
         if (onInlet_(globalPos))
             values[conti0EqIdx] = -0.04; // kg/(m*s)
         return values;
diff --git a/test/porousmediumflow/richards/implicit/richardsniconductionproblem.hh b/test/porousmediumflow/richards/implicit/richardsniconductionproblem.hh
index 4f1510a00a7738b2393547d00ae97f42335923af..a2591ef4927423ff5d411f5a7249cd0d0c239474 100644
--- a/test/porousmediumflow/richards/implicit/richardsniconductionproblem.hh
+++ b/test/porousmediumflow/richards/implicit/richardsniconductionproblem.hh
@@ -101,6 +101,7 @@ class RichardsNIConductionProblem :public PorousMediumFlowProblem<TypeTag>
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);
     using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
     using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
@@ -256,9 +257,9 @@ public:
      * in normal direction of each component. Negative values mean
      * influx.
      */
-    PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
+    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
     {
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
         return values;
     }
 
diff --git a/test/porousmediumflow/richards/implicit/richardsniconvectionproblem.hh b/test/porousmediumflow/richards/implicit/richardsniconvectionproblem.hh
index ed504f860435111802c2fd8b4089e591541e5d18..cfeedd06f892b3c8b9877903738f3e1f0ba66aa2 100644
--- a/test/porousmediumflow/richards/implicit/richardsniconvectionproblem.hh
+++ b/test/porousmediumflow/richards/implicit/richardsniconvectionproblem.hh
@@ -105,6 +105,7 @@ class RichardsNIConvectionProblem : public PorousMediumFlowProblem<TypeTag>
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
     using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);
@@ -266,12 +267,12 @@ public:
      * \param scvf The subcontrolvolume face
      *  Negative values mean influx.
      */
-    PrimaryVariables neumann(const Element &element,
+    NumEqVector neumann(const Element &element,
                              const FVElementGeometry& fvGeometry,
                              const ElementVolumeVariables& elemVolVars,
                              const SubControlVolumeFace& scvf) const
     {
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
         const auto globalPos = scvf.ipGlobal();
 
         if(globalPos[0] < eps_)
diff --git a/test/porousmediumflow/richardsnc/implicit/richardswelltracerproblem.hh b/test/porousmediumflow/richardsnc/implicit/richardswelltracerproblem.hh
index 496f2af5d5b40a67a8a7f28af60a995141a8d600..6fafc3690b87527766f022c333f328c901eee19b 100644
--- a/test/porousmediumflow/richardsnc/implicit/richardswelltracerproblem.hh
+++ b/test/porousmediumflow/richardsnc/implicit/richardswelltracerproblem.hh
@@ -101,6 +101,7 @@ class RichardsWellTracerProblem : public PorousMediumFlowProblem<TypeTag>
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using PointSource = typename GET_PROP_TYPE(TypeTag, PointSource);
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
@@ -250,8 +251,8 @@ public:
      *
      * \param globalPos The position for which the Neumann value is set
      */
-    PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
-    { return PrimaryVariables(0.0); }
+    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
+    { return NumEqVector(0.0); }
 
     /*!
      * \name Volume terms
diff --git a/test/porousmediumflow/tracer/multicomp/maxwellstefantestproblem.hh b/test/porousmediumflow/tracer/multicomp/maxwellstefantestproblem.hh
index 8b97ab35673e50749b41516fd4fee52236aab094..a4398626fe2edd130c13d2f6b3c6bdcb54ad360d 100644
--- a/test/porousmediumflow/tracer/multicomp/maxwellstefantestproblem.hh
+++ b/test/porousmediumflow/tracer/multicomp/maxwellstefantestproblem.hh
@@ -181,6 +181,7 @@ class MaxwellStefanTestProblem : public PorousMediumFlowProblem<TypeTag>
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
     using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
@@ -339,8 +340,8 @@ public:
      * \param globalPos The position for which the bc type should be evaluated
      * The units must be according to either using mole or mass fractions. (mole/(m^2*s) or kg/(m^2*s))
      */
-    PrimaryVariables neumannAtPos(const GlobalPosition& globalPos) const
-    { return PrimaryVariables(0.0); }
+    NumEqVector neumannAtPos(const GlobalPosition& globalPos) const
+    { return NumEqVector(0.0); }
 
     // \}