diff --git a/exercises/exercise-basic/injection2p2cproblem.hh b/exercises/exercise-basic/injection2p2cproblem.hh
index f564bffe04205ab4d8a137a1fd7a255c22bf9886..e3ba1167b5bec00b52f9fcfa2b2fdebac3fa7a31 100644
--- a/exercises/exercise-basic/injection2p2cproblem.hh
+++ b/exercises/exercise-basic/injection2p2cproblem.hh
@@ -197,10 +197,10 @@ public:
      *
      * \param globalPos The position of the integration point of the boundary segment.
      */
-    PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
+    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
     {
         // initialize values to zero, i.e. no-flow Neumann boundary conditions
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
 
         // if we are inside the injection zone set inflow Neumann boundary conditions
         if (time_ < injectionDuration_
diff --git a/exercises/exercise-basic/injection2pniproblem.hh b/exercises/exercise-basic/injection2pniproblem.hh
index 7642c96a53db03c2eb80ba88f526f50fed8890ed..3b90fd9f8a3de38d342c6d3290ea83af91a8e765 100644
--- a/exercises/exercise-basic/injection2pniproblem.hh
+++ b/exercises/exercise-basic/injection2pniproblem.hh
@@ -195,10 +195,10 @@ public:
      *
      * \param globalPos The position of the integration point of the boundary segment.
      */
-    PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
+    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
     {
         // initialize values to zero, i.e. no-flow Neumann boundary conditions
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
 
         // if we are inside the injection zone set inflow Neumann boundary conditions
         if (time_ < injectionDuration_
diff --git a/exercises/exercise-basic/injection2pproblem.hh b/exercises/exercise-basic/injection2pproblem.hh
index a4f762b51cb01f347ca5cb2dd2815dcbcc2cae78..5dc9a4b6f969a49c96fd0c697988eb17e3243c7e 100644
--- a/exercises/exercise-basic/injection2pproblem.hh
+++ b/exercises/exercise-basic/injection2pproblem.hh
@@ -197,10 +197,10 @@ public:
      *
      * \param globalPos The position of the integration point of the boundary segment.
      */
-    PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
+    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
     {
         // initialize values to zero, i.e. no-flow Neumann boundary conditions
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
 
         // if we are inside the injection zone set inflow Neumann boundary conditions
         // using < boundary + eps_ or > boundary - eps_ is safer for floating point comparisons
diff --git a/exercises/exercise-fluidsystem/2p2cproblem.hh b/exercises/exercise-fluidsystem/2p2cproblem.hh
index 0c3082ca1650f529e3a8f3c796f7325208f36bd8..e8aef797d6dd3d4cdccc664561cf2826d44868bb 100644
--- a/exercises/exercise-fluidsystem/2p2cproblem.hh
+++ b/exercises/exercise-fluidsystem/2p2cproblem.hh
@@ -110,6 +110,7 @@ class ExerciseFluidsystemProblemTwoPTwoC : public PorousMediumFlowProblem<TypeTa
     using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
     using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
 
 public:
     ExerciseFluidsystemProblemTwoPTwoC(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
@@ -176,10 +177,10 @@ public:
      *
      * \param globalPos The position of the integration point of the boundary segment.
      */
-    PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
+    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
     {
         // initialize values to zero, i.e. no-flow Neumann boundary conditions
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
 
         Scalar up = this->fvGridGeometry().bBoxMax()[dimWorld-1];
         // extraction of oil (30 g/m/s) on a segment of the upper boundary
@@ -236,10 +237,10 @@ public:
      * \brief Returns the source term
      * \param globalPos The global position
      */
-    PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
     {
         // we do not define any sources
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
         return values;
     }
 
diff --git a/exercises/exercise-fluidsystem/2pproblem.hh b/exercises/exercise-fluidsystem/2pproblem.hh
index 1acfb0ff916b27b4089c92608c6676ebfed3d56a..4462997fa1ed2095e1d45ee9966d0d25e8967cf1 100644
--- a/exercises/exercise-fluidsystem/2pproblem.hh
+++ b/exercises/exercise-fluidsystem/2pproblem.hh
@@ -133,6 +133,7 @@ class ExerciseFluidsystemProblemTwoP : public PorousMediumFlowProblem<TypeTag>
     using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using FluidState = GetPropType<TypeTag, Properties::FluidState>;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
 
     enum {
         waterPressureIdx = Indices::pressureIdx,
@@ -218,10 +219,10 @@ 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
     {
         // initialize values to zero, i.e. no-flow Neumann boundary conditions
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
 
         Scalar up = this->fvGridGeometry().bBoxMax()[dimWorld-1];
 
@@ -274,10 +275,10 @@ public:
      *
      * \param globalPos The global position
      */
-    PrimaryVariables sourceAtPos(const GlobalPosition& globalPos) const
+    NumEqVector sourceAtPos(const GlobalPosition& globalPos) const
     {
         // we define no source term here
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
         return values;
     }
 
diff --git a/exercises/exercise-grids/injection2pproblem.hh b/exercises/exercise-grids/injection2pproblem.hh
index 910b2ff8a5dfebe6ca14442bd7ae7ddc5f3190f6..af036114eeab1199a567d0ba16774a76fe6b3a78 100644
--- a/exercises/exercise-grids/injection2pproblem.hh
+++ b/exercises/exercise-grids/injection2pproblem.hh
@@ -108,6 +108,7 @@ class InjectionProblem2P : public PorousMediumFlowProblem<TypeTag>
     using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
     using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
 
     enum { dimWorld = GridView::dimensionworld };
     using Element = typename GridView::template Codim<0>::Entity;
@@ -199,10 +200,10 @@ public:
      *
      * \param globalPos The position of the integration point of the boundary segment.
      */
-    PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
+    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
     {
         // initialize values to zero, i.e. no-flow Neumann boundary conditions
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
 
         // if we are inside the injection zone set inflow Neumann boundary conditions
         // using < boundary + eps_ or > boundary - eps_ is safer for floating point comparisons
diff --git a/exercises/exercise-runtimeparams/injection2pproblem.hh b/exercises/exercise-runtimeparams/injection2pproblem.hh
index ea353268dfbe27f1bdffed38d1b0d2985dd416d0..6f6500e9e10648c5ecd7f4cf3215a8113ceb1808 100644
--- a/exercises/exercise-runtimeparams/injection2pproblem.hh
+++ b/exercises/exercise-runtimeparams/injection2pproblem.hh
@@ -105,6 +105,7 @@ class InjectionProblem2P : public PorousMediumFlowProblem<TypeTag>
     using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
     using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
 
     enum { dimWorld = GridView::dimensionworld };
     using Element = typename GridView::template Codim<0>::Entity;
@@ -199,10 +200,10 @@ public:
      *
      * \param globalPos The position of the integration point of the boundary segment.
      */
-    PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
+    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
     {
         // initialize values to zero, i.e. no-flow Neumann boundary conditions
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
 
         // if we are inside the injection zone set inflow Neumann boundary conditions
         // using < boundary + eps_ or > boundary - eps_ is safer for floating point comparisons
diff --git a/exercises/solution/exercise-basic/injection2pniproblem.hh b/exercises/solution/exercise-basic/injection2pniproblem.hh
index 8d971152777ecdd89e5b7d1b531bb6578e5e72c2..e722c06df1a10dab7cf5ca2174b003522dc0bba1 100644
--- a/exercises/solution/exercise-basic/injection2pniproblem.hh
+++ b/exercises/solution/exercise-basic/injection2pniproblem.hh
@@ -190,10 +190,10 @@ public:
      *
      * \param globalPos The position of the integration point of the boundary segment.
      */
-    PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
+    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
     {
         // initialize values to zero, i.e. no-flow Neumann boundary conditions
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
 
         // if we are inside the injection zone set inflow Neumann boundary conditions
         if (time_ < injectionDuration_
diff --git a/exercises/solution/exercise-fluidsystem/2p2cproblem.hh b/exercises/solution/exercise-fluidsystem/2p2cproblem.hh
index 25ae5d6c9939ea9e2509f2ecf07d7455944a8ce1..b015f84f8184be1acd0cd9e0fb0d24701a9faa6b 100644
--- a/exercises/solution/exercise-fluidsystem/2p2cproblem.hh
+++ b/exercises/solution/exercise-fluidsystem/2p2cproblem.hh
@@ -110,6 +110,7 @@ class ExerciseFluidsystemProblemTwoPTwoC : public PorousMediumFlowProblem<TypeTa
     using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
     using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
 
 public:
     ExerciseFluidsystemProblemTwoPTwoC(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
@@ -181,10 +182,10 @@ 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
     {
         // initialize values to zero, i.e. no-flow Neumann boundary conditions
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
 
         Scalar up = this->fvGridGeometry().bBoxMax()[dimWorld-1];
         // extraction of oil (30 g/m/s) on a segment of the upper boundary
@@ -245,10 +246,10 @@ 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
     {
         // we do not define any sources
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
         return values;
     }
 
diff --git a/exercises/solution/exercise-fluidsystem/2pproblem.hh b/exercises/solution/exercise-fluidsystem/2pproblem.hh
index fa461792f0ae45a1262d4add2c6c52bcc287db20..92726d10888620f36fcf940709b3d29dfbfababe 100644
--- a/exercises/solution/exercise-fluidsystem/2pproblem.hh
+++ b/exercises/solution/exercise-fluidsystem/2pproblem.hh
@@ -133,6 +133,7 @@ class ExerciseFluidsystemProblemTwoP : public PorousMediumFlowProblem<TypeTag>
     using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using FluidState = GetPropType<TypeTag, Properties::FluidState>;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
 
     enum {
         waterPressureIdx = Indices::pressureIdx,
@@ -221,10 +222,10 @@ 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
     {
         // initialize values to zero, i.e. no-flow Neumann boundary conditions
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
 
         Scalar up = this->fvGridGeometry().bBoxMax()[dimWorld-1];
 
@@ -279,10 +280,10 @@ 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
     {
         // we define no source term here
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
         return values;
     }
 
diff --git a/exercises/solution/exercise-grids/injection2pproblem.hh b/exercises/solution/exercise-grids/injection2pproblem.hh
index 101c0379150b600d40dbe7404ec669f03e0b3da2..df0e9d08836912652fcd8c7aaa105188e7546164 100644
--- a/exercises/solution/exercise-grids/injection2pproblem.hh
+++ b/exercises/solution/exercise-grids/injection2pproblem.hh
@@ -116,6 +116,7 @@ class InjectionProblem2P : public PorousMediumFlowProblem<TypeTag>
     using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
     using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
 
     enum { dimWorld = GridView::dimensionworld };
     using Element = typename GridView::template Codim<0>::Entity;
@@ -207,10 +208,10 @@ public:
      *
      * \param globalPos The position of the integration point of the boundary segment.
      */
-    PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
+    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
     {
         // initialize values to zero, i.e. no-flow Neumann boundary conditions
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
 
         // if we are inside the injection zone set inflow Neumann boundary conditions
         // using < boundary + eps_ or > boundary - eps_ is safer for floating point comparisons
diff --git a/exercises/solution/exercise-runtimeparams/injection2pproblem.hh b/exercises/solution/exercise-runtimeparams/injection2pproblem.hh
index b2d79304a3f0a2778f1b075ad5b7b41241c7d30f..b4ce84ef654e9af2889b829c720c146dcce259da 100644
--- a/exercises/solution/exercise-runtimeparams/injection2pproblem.hh
+++ b/exercises/solution/exercise-runtimeparams/injection2pproblem.hh
@@ -105,6 +105,7 @@ class InjectionProblem2P : public PorousMediumFlowProblem<TypeTag>
     using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
     using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
 
     enum { dimWorld = GridView::dimensionworld };
     using Element = typename GridView::template Codim<0>::Entity;
@@ -205,10 +206,10 @@ public:
      *
      * \param globalPos The position of the integration point of the boundary segment.
      */
-    PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
+    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
     {
         // initialize values to zero, i.e. no-flow Neumann boundary conditions
-        PrimaryVariables values(0.0);
+        NumEqVector values(0.0);
 
         // if we are inside the injection zone set inflow Neumann boundary conditions
         // using < boundary + eps_ or > boundary - eps_ is safer for floating point comparisons