diff --git a/dumux/assembly/fvassembler.hh b/dumux/assembly/fvassembler.hh
index d26121d52cba464655d356239284f1bdf0405378..dddfe595b6da961d84af08cff945f6250cfd8f10 100644
--- a/dumux/assembly/fvassembler.hh
+++ b/dumux/assembly/fvassembler.hh
@@ -51,7 +51,6 @@ template<class TypeTag, DiffMethod diffMethod, bool isImplicit = true>
 class FVAssembler
 {
     using GridView = GetPropType<TypeTag, Properties::GridView>;
-    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
     using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>;
     using Element = typename GridView::template Codim<0>::Entity;
     using TimeLoop = TimeLoopBase<GetPropType<TypeTag, Properties::Scalar>>;
@@ -69,6 +68,7 @@ public:
     using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
     using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
     using Problem = GetPropType<TypeTag, Properties::Problem>;
+    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
 
     using ResidualType = SolutionVector;
 
diff --git a/dumux/assembly/staggeredfvassembler.hh b/dumux/assembly/staggeredfvassembler.hh
index 9e8948ac0930998829d061b4da0a97da3f01bf50..c5f17e9257381521ba00d87d82fa95f2e04c1c91 100644
--- a/dumux/assembly/staggeredfvassembler.hh
+++ b/dumux/assembly/staggeredfvassembler.hh
@@ -61,11 +61,11 @@ class StaggeredFVAssembler: public MultiDomainFVAssembler<StaggeredMultiDomainTr
                                               diffMethod>;
 
     using Problem = GetPropType<TypeTag, Properties::Problem>;
-    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
     using TimeLoop = TimeLoopBase<GetPropType<TypeTag, Properties::Scalar>>;
 
 public:
     using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
     using CouplingManager = typename ParentType::CouplingManager;
 
     //! The constructor for stationary problems
diff --git a/dumux/common/properties.hh b/dumux/common/properties.hh
index acde51138f8bcf8c06c658337a3c8bcb7fa38fa2..7d3e84f123713815b30092e8eba5888aed7e1511 100644
--- a/dumux/common/properties.hh
+++ b/dumux/common/properties.hh
@@ -161,8 +161,6 @@ struct SolidSystem { using type = UndefinedProperty; };
 template<class TypeTag, class MyTypeTag>
 struct SolidState { using type = UndefinedProperty; };                           //!< The type of the solid state to use
 template<class TypeTag, class MyTypeTag>
-struct PrimaryVariableSwitch { using type = UndefinedProperty; };               //!< The primary variable switch needed for compositional models
-template<class TypeTag, class MyTypeTag>
 struct EffectiveDiffusivityModel { using type = UndefinedProperty; };           //!< The employed model for the computation of the effective diffusivity
 template<class TypeTag, class MyTypeTag>
 struct ThermalConductivityModel { using type = UndefinedProperty; };            //!< Model to be used for the calculation of the effective conductivity
diff --git a/dumux/multidomain/CMakeLists.txt b/dumux/multidomain/CMakeLists.txt
index 28588a83ff2513e73cf3869f4abb65ba5e3159d5..c983c6a65566e81c51d0770d7f099b4a2a4d4845 100644
--- a/dumux/multidomain/CMakeLists.txt
+++ b/dumux/multidomain/CMakeLists.txt
@@ -7,7 +7,6 @@ couplingjacobianpattern.hh
 couplingmanager.hh
 fvassembler.hh
 newtonsolver.hh
-privarswitchnewtonsolver.hh
 staggeredcouplingmanager.hh
 staggeredtraits.hh
 subdomainboxlocalassembler.hh
diff --git a/dumux/multidomain/fvassembler.hh b/dumux/multidomain/fvassembler.hh
index 34c2e4a9ca72bb013c654b816758c293a8dd6f73..13e80fa458b1e8edd00e2b90c5d639096670954b 100644
--- a/dumux/multidomain/fvassembler.hh
+++ b/dumux/multidomain/fvassembler.hh
@@ -27,6 +27,7 @@
 #define DUMUX_MULTIDOMAIN_FV_ASSEMBLER_HH
 
 #include <type_traits>
+#include <tuple>
 
 #include <dune/common/hybridutilities.hh>
 #include <dune/istl/matrixindexset.hh>
@@ -67,6 +68,9 @@ public:
     template<std::size_t id>
     using LocalResidual = GetPropType<SubDomainTypeTag<id>, Properties::LocalResidual>;
 
+    template<std::size_t id>
+    using GridVariables = typename std::tuple_element_t<id, typename MDTraits::GridVariablesTuple>::element_type;
+
     using JacobianMatrix = typename MDTraits::JacobianMatrix;
     using SolutionVector = typename MDTraits::SolutionVector;
     using ResidualType = SolutionVector;
@@ -337,12 +341,12 @@ public:
 
     //! the grid variables of domain i
     template<std::size_t i>
-    auto& gridVariables(Dune::index_constant<i> domainId)
+    GridVariables<i>& gridVariables(Dune::index_constant<i> domainId)
     { return *std::get<domainId>(gridVariablesTuple_); }
 
     //! the grid variables of domain i
     template<std::size_t i>
-    const auto& gridVariables(Dune::index_constant<i> domainId) const
+    const GridVariables<i>& gridVariables(Dune::index_constant<i> domainId) const
     { return *std::get<domainId>(gridVariablesTuple_); }
 
     //! the coupling manager
diff --git a/dumux/multidomain/newtonsolver.hh b/dumux/multidomain/newtonsolver.hh
index 8a661a8475c4dc12e7baee4a1fe25190839bb23e..e63fb4fd77de82d57a772b13efeed42b53e6bc12 100644
--- a/dumux/multidomain/newtonsolver.hh
+++ b/dumux/multidomain/newtonsolver.hh
@@ -29,6 +29,15 @@
 #include <dumux/nonlinear/newtonsolver.hh>
 
 namespace Dumux {
+namespace Detail {
+
+template<class Assembler, class Index>
+using DetectPVSwitchMultiDomain = typename Assembler::template GridVariables<Index::value>::VolumeVariables::PrimaryVariableSwitch;
+
+template<class Assembler, std::size_t i>
+using GetPVSwitchMultiDomain = Dune::Std::detected_or<int, DetectPVSwitchMultiDomain, Assembler, Dune::index_constant<i>>;
+
+} // end namespace Detail
 
 /*!
  * \ingroup Nonlinear
@@ -43,6 +52,16 @@ class MultiDomainNewtonSolver: public NewtonSolver<Assembler, LinearSolver, Reas
     using ParentType = NewtonSolver<Assembler, LinearSolver, Reassembler, Comm>;
     using SolutionVector = typename Assembler::ResidualType;
 
+    template<std::size_t i>
+    using PrimaryVariableSwitch = typename Detail::GetPVSwitchMultiDomain<Assembler, i>::type;
+
+    template<std::size_t i>
+    using HasPriVarsSwitch = typename Detail::GetPVSwitchMultiDomain<Assembler, i>::value_t; // std::true_type or std::false_type
+
+    template<std::size_t i>
+    using PrivarSwitchPtr = std::unique_ptr<PrimaryVariableSwitch<i>>;
+    using PriVarSwitchPtrTuple = typename Assembler::Traits::template MakeTuple<PrivarSwitchPtr>;
+
 public:
 
     /*!
@@ -55,7 +74,17 @@ public:
                             const std::string& paramGroup = "")
     : ParentType(assembler, linearSolver, comm, paramGroup)
     , couplingManager_(couplingManager)
-    {}
+    {
+        using namespace Dune::Hybrid;
+        forEach(std::make_index_sequence<Assembler::Traits::numSubDomains>{}, [&](auto&& id)
+        {
+            const int priVarSwitchVerbosity = getParamFromGroup<int>(paramGroup, "PrimaryVariableSwitch.Verbosity", 1);
+            using PVSwitch = PrimaryVariableSwitch<std::decay_t<decltype(id)>::value>;
+            elementAt(priVarSwitches_, id) = std::make_unique<PVSwitch>(priVarSwitchVerbosity);
+        });
+
+        priVarsSwitchedInLastIteration_.fill(false);
+    }
 
     /*!
      * \brief Indicates the beginning of a Newton iteration.
@@ -66,6 +95,36 @@ public:
         couplingManager_->updateSolution(uCurrentIter);
     }
 
+    /*!
+     *
+     * \brief Called before the Newton method is applied to an
+     *        non-linear system of equations.
+     *
+     * \param u The initial solution
+     */
+    void newtonBegin(const SolutionVector &u) override
+    {
+        ParentType::newtonBegin(u);
+
+        using namespace Dune::Hybrid;
+        forEach(std::make_index_sequence<Assembler::Traits::numSubDomains>{}, [&](auto&& id)
+        {
+            resetPriVarSwitch_(u[id].size(), id, HasPriVarsSwitch<id>{});
+        });
+    }
+
+    /*!
+     * \brief Returns true if the error of the solution is below the
+     *        tolerance.
+     */
+    bool newtonConverged() const override
+    {
+        if (Dune::any_true(priVarsSwitchedInLastIteration_))
+            return false;
+
+        return ParentType::newtonConverged();
+    }
+
 
     /*!
      * \brief Indicates that one Newton iteration was finished.
@@ -76,12 +135,80 @@ public:
      */
     void newtonEndStep(SolutionVector &uCurrentIter, const SolutionVector &uLastIter) override
     {
+        using namespace Dune::Hybrid;
+        forEach(std::make_index_sequence<Assembler::Traits::numSubDomains>{}, [&](auto&& id)
+        {
+            invokePriVarSwitch_(uCurrentIter[id], id, HasPriVarsSwitch<id>{});
+        });
+
         ParentType::newtonEndStep(uCurrentIter, uLastIter);
         couplingManager_->updateSolution(uCurrentIter);
     }
 
 private:
+
+    /*!
+     * \brief Reset the privar switch state, noop if there is no priVarSwitch
+     */
+    template<std::size_t i>
+    void resetPriVarSwitch_(const std::size_t numDofs, Dune::index_constant<i> id, std::false_type) {}
+
+    /*!
+     * \brief Switch primary variables if necessary
+     */
+    template<std::size_t i>
+    void resetPriVarSwitch_(const std::size_t numDofs, Dune::index_constant<i> id, std::true_type)
+    {
+        using namespace Dune::Hybrid;
+        elementAt(priVarSwitches_, id)->reset(numDofs);
+        priVarsSwitchedInLastIteration_[i] = false;
+    }
+
+    /*!
+     * \brief Switch primary variables if necessary, noop if there is no priVarSwitch
+     */
+    template<class SubSol, std::size_t i>
+    void invokePriVarSwitch_(SubSol&, Dune::index_constant<i> id, std::false_type) {}
+
+    /*!
+     * \brief Switch primary variables if necessary
+     */
+    template<class SubSol, std::size_t i>
+    void invokePriVarSwitch_(SubSol& uCurrentIter, Dune::index_constant<i> id, std::true_type)
+    {
+        // update the variable switch (returns true if the pri vars at at least one dof were switched)
+        // for disabled grid variable caching
+        const auto& fvGridGeometry = this->assembler().fvGridGeometry(id);
+        const auto& problem = this->assembler().problem(id);
+        auto& gridVariables = this->assembler().gridVariables(id);
+
+        using namespace Dune::Hybrid;
+        auto& priVarSwitch = *elementAt(priVarSwitches_, id);
+
+        // invoke the primary variable switch
+        priVarsSwitchedInLastIteration_[i] = priVarSwitch.update(uCurrentIter, gridVariables,
+                                                                 problem, fvGridGeometry);
+
+        if (priVarsSwitchedInLastIteration_[i])
+        {
+            for (const auto& element : elements(fvGridGeometry.gridView()))
+            {
+                // if the volume variables are cached globally, we need to update those where the primary variables have been switched
+                priVarSwitch.updateSwitchedVolVars(problem, element, fvGridGeometry, gridVariables, uCurrentIter);
+
+                // if the flux variables are cached globally, we need to update those where the primary variables have been switched
+                priVarSwitch.updateSwitchedFluxVarsCache(problem, element, fvGridGeometry, gridVariables, uCurrentIter);
+            }
+        }
+    }
+
+    //! the coupling manager
     std::shared_ptr<CouplingManager> couplingManager_;
+
+    //! the class handling the primary variable switch
+    PriVarSwitchPtrTuple priVarSwitches_;
+    //! if we switched primary variables in the last iteration
+    std::array<bool, Assembler::Traits::numSubDomains> priVarsSwitchedInLastIteration_;
 };
 
 } // end namespace Dumux
diff --git a/dumux/multidomain/privarswitchnewtonsolver.hh b/dumux/multidomain/privarswitchnewtonsolver.hh
deleted file mode 100644
index f4961e7fa85ad7c60c343c00b65611585ebfca28..0000000000000000000000000000000000000000
--- a/dumux/multidomain/privarswitchnewtonsolver.hh
+++ /dev/null
@@ -1,160 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/****************************************************************************
- *   See the file COPYING for full copying permissions.                      *
- *                                                                           *
- *   This program is free software: you can redistribute it and/or modify    *
- *   it under the terms of the GNU General Public License as published by    *
- *   the Free Software Foundation, either version 2 of the License, or       *
- *   (at your option) any later version.                                     *
- *                                                                           *
- *   This program is distributed in the hope that it will be useful,         *
- *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
- *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
- *   GNU General Public License for more details.                            *
- *                                                                           *
- *   You should have received a copy of the GNU General Public License       *
- *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
- *****************************************************************************/
-/*!
- * \file
- * \ingroup Nonlinear
- * \ingroup MultiDomain
- * \copydoc Dumux::MultiDomainPriVarSwitchNewtonSolver
- */
-#ifndef DUMUX_MULTIDOMAIN_PRIVARSWITCH_NEWTON_SOLVER_HH
-#define DUMUX_MULTIDOMAIN_PRIVARSWITCH_NEWTON_SOLVER_HH
-
-#include <memory>
-
-#include <dumux/common/typetraits/utility.hh>
-#include <dumux/discretization/method.hh>
-#include <dumux/discretization/elementsolution.hh>
-#include <dumux/multidomain/newtonsolver.hh>
-
-namespace Dumux {
-
-/*!
- * \ingroup Nonlinear
- * \ingroup MultiDomain
- * \brief A newton solver that handles primary variable switches for multi domain problems
- */
-template <class Assembler, class LinearSolver, class CouplingManager, class PrivarSwitchTypeTuple,
-          class Reassembler = DefaultPartialReassembler,
-          class Comm = Dune::CollectiveCommunication<Dune::MPIHelper::MPICommunicator> >
-class MultiDomainPriVarSwitchNewtonSolver : public MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager, Reassembler, Comm>
-{
-    using Scalar =  typename Assembler::Scalar;
-    using ParentType = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager, Reassembler, Comm>;
-    using SolutionVector = typename Assembler::ResidualType;
-
-    template<std::size_t i> using PriVarSwitch = std::tuple_element_t<i, PrivarSwitchTypeTuple>;
-    template<std::size_t i> using PrivarSwitchPtr = std::unique_ptr<PriVarSwitch<i>>;
-    using PriVarSwitchPtrTuple = typename Assembler::Traits::template MakeTuple<PrivarSwitchPtr>;
-
-public:
-    using ParentType::ParentType;
-
-    /*!
-     * \brief Returns true if the error of the solution is below the
-     *        tolerance.
-     */
-    bool newtonConverged() const override
-    {
-        if (std::any_of(switchedInLastIteration_.begin(), switchedInLastIteration_.end(), [](bool i){ return i;}))
-            return false;
-
-        return ParentType::newtonConverged();
-    }
-
-    /*!
-     *
-     * \brief Called before the Newton method is applied to an
-     *        non-linear system of equations.
-     *
-     * \param u The initial solution
-     */
-    void newtonBegin(const SolutionVector &u) override
-    {
-        ParentType::newtonBegin(u);
-        switchedInLastIteration_.fill(false);
-
-        using namespace Dune::Hybrid;
-        forEach(std::make_index_sequence<Assembler::Traits::numSubDomains>{}, [&](auto&& id)
-        {
-            using PVSwitch = PriVarSwitch<std::decay_t<decltype(id)>::value>;
-            elementAt(priVarSwitches_, id) = std::make_unique<PVSwitch>(u[id].size());
-        });
-    }
-
-    /*!
-     * \brief Indicates that one Newton iteration was finished.
-     *
-     * \param assembler The jacobian assembler
-     * \param uCurrentIter The solution after the current Newton iteration
-     * \param uLastIter The solution at the beginning of the current Newton iteration
-     */
-    void newtonEndStep(SolutionVector &uCurrentIter,
-                       const SolutionVector &uLastIter) override
-    {
-        ParentType::newtonEndStep(uCurrentIter, uLastIter);
-
-        // update the variable switch (returns true if the pri vars at at least one dof were switched)
-        // for disabled grid variable caching
-        auto& assembler = this->assembler();
-
-        using namespace Dune::Hybrid;
-        forEach(std::make_index_sequence<Assembler::Traits::numSubDomains>{}, [&](auto&& id)
-        {
-            const auto& fvGridGeometry = assembler.fvGridGeometry(id);
-            const auto& problem = assembler.problem(id);
-            auto& gridVariables = assembler.gridVariables(id);
-            auto& priVarSwitch = elementAt(priVarSwitches_, id);
-
-            // invoke the primary variable switch
-            switchedInLastIteration_[id] = priVarSwitch->update(uCurrentIter[id], gridVariables,
-                                                                               problem, fvGridGeometry);
-
-            if(switchedInLastIteration_[id])
-            {
-                for (const auto& element : elements(fvGridGeometry.gridView()))
-                {
-                    // if the volume variables are cached globally, we need to update those where the primary variables have been switched
-                    priVarSwitch->updateSwitchedVolVars(problem, element, fvGridGeometry, gridVariables, uCurrentIter);
-
-                    // if the flux variables are cached globally, we need to update those where the primary variables have been switched
-                    priVarSwitch->updateSwitchedFluxVarsCache(problem, element, fvGridGeometry, gridVariables, uCurrentIter);
-                }
-            }
-        });
-    }
-
-    /*!
-     * \brief Called if the Newton method ended
-     *        (not known yet if we failed or succeeded)
-     */
-    void newtonEnd() override
-    {
-        ParentType::newtonEnd();
-
-        // in any way, we have to reset the switch flag
-        switchedInLastIteration_.fill(false);
-
-        // free some memory
-        using namespace Dune::Hybrid;
-        forEach(std::make_index_sequence<Assembler::Traits::numSubDomains>{}, [&](auto&& id)
-        {
-            elementAt(priVarSwitches_, id).release();
-        });
-    }
-
-private:
-    //! the class handling the primary variable switch
-    PriVarSwitchPtrTuple priVarSwitches_;
-    //! if we switched primary variables in the last iteration
-    std::array<bool, Assembler::Traits::numSubDomains> switchedInLastIteration_;
-};
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/nonlinear/CMakeLists.txt b/dumux/nonlinear/CMakeLists.txt
index ce718221685cb34e3c82a14a0b06b71f77481a39..1ab033a40e7803299a74a595dde95bda4ee19118 100644
--- a/dumux/nonlinear/CMakeLists.txt
+++ b/dumux/nonlinear/CMakeLists.txt
@@ -1,5 +1,4 @@
 install(FILES
 newtonconvergencewriter.hh
 newtonsolver.hh
-privarswitchnewtonsolver.hh
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/nonlinear)
diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh
index 9df3c0bfd4687dce944570c89323301d915eee61..50b124211369dafe9122865172438832ee2d3f73 100644
--- a/dumux/nonlinear/newtonsolver.hh
+++ b/dumux/nonlinear/newtonsolver.hh
@@ -29,11 +29,13 @@
 #include <cmath>
 #include <memory>
 #include <iostream>
+#include <type_traits>
 
 #include <dune/common/timer.hh>
 #include <dune/common/exceptions.hh>
 #include <dune/common/parallel/mpicollectivecommunication.hh>
 #include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/std/type_traits.hh>
 #include <dune/istl/bvector.hh>
 #include <dune/istl/multitypeblockvector.hh>
 
@@ -48,6 +50,7 @@
 #include "newtonconvergencewriter.hh"
 
 namespace Dumux {
+namespace Detail {
 
 //! helper struct detecting if an assembler supports partial reassembly
 struct supportsPartialReassembly
@@ -59,6 +62,15 @@ struct supportsPartialReassembly
     {}
 };
 
+//! helper aliases to extract a primary variable switch from the VolumeVariables (if defined, yields int otherwise)
+template<class Assembler>
+using DetectPVSwitch = typename Assembler::GridVariables::VolumeVariables::PrimaryVariableSwitch;
+
+template<class Assembler>
+using GetPVSwitch = Dune::Std::detected_or<int, DetectPVSwitch, Assembler>;
+
+} // end namespace Detail
+
 /*!
  * \ingroup Nonlinear
  * \brief An implementation of a Newton solver
@@ -79,6 +91,10 @@ class NewtonSolver
     using SolutionVector = typename Assembler::ResidualType;
     using ConvergenceWriter = ConvergenceWriterInterface<SolutionVector>;
 
+    using PrimaryVariableSwitch = typename Detail::GetPVSwitch<Assembler>::type;
+    using HasPriVarsSwitch = typename Detail::GetPVSwitch<Assembler>::value_t; // std::true_type or std::false_type
+    static constexpr bool hasPriVarsSwitch() { return HasPriVarsSwitch::value; };
+
 public:
 
     using Communication = Comm;
@@ -108,6 +124,12 @@ public:
         // initialize the partial reassembler
         if (enablePartialReassembly_)
             partialReassembler_ = std::make_unique<Reassembler>(*assembler_);
+
+        if (hasPriVarsSwitch())
+        {
+            const int priVarSwitchVerbosity = getParamFromGroup<int>(paramGroup, "PrimaryVariableSwitch.Verbosity", 1);
+            priVarSwitch_ = std::make_unique<PrimaryVariableSwitch>(priVarSwitchVerbosity);
+        }
     }
 
     virtual ~NewtonSolver() {}
@@ -233,6 +255,7 @@ public:
     virtual void newtonBegin(const SolutionVector& u)
     {
         numSteps_ = 0;
+        resetPriVarSwitch_(u.size(), HasPriVarsSwitch{});
     }
 
     /*!
@@ -445,6 +468,8 @@ public:
     virtual void newtonEndStep(SolutionVector &uCurrentIter,
                                const SolutionVector &uLastIter)
     {
+        invokePriVarSwitch_(uCurrentIter, HasPriVarsSwitch{});
+
         ++numSteps_;
 
         if (verbose_)
@@ -487,6 +512,11 @@ public:
      */
     virtual bool newtonConverged() const
     {
+        // in case the model has a priVar switch and some some primary variables
+        // actually switched their state in the last iteration, enforce another iteration
+        if (priVarsSwitchedInLastIteration_)
+            return false;
+
         if (enableShiftCriterion_ && !enableResidualCriterion_)
         {
             return shift_ <= shiftTolerance_;
@@ -621,6 +651,53 @@ protected:
     Assembler& assembler()
     { return *assembler_; }
 
+    /*!
+     * \brief Reset the privar switch state, noop if there is no priVarSwitch
+     */
+    void resetPriVarSwitch_(const std::size_t numDofs, std::false_type) {}
+
+    /*!
+     * \brief Reset the privar switch state
+     */
+    void resetPriVarSwitch_(const std::size_t numDofs, std::true_type)
+    {
+        priVarSwitch_->reset(numDofs);
+        priVarsSwitchedInLastIteration_ = false;
+    }
+
+    /*!
+     * \brief Switch primary variables if necessary, noop if there is no priVarSwitch
+     */
+    void invokePriVarSwitch_(SolutionVector&, std::false_type) {}
+
+    /*!
+     * \brief Switch primary variables if necessary
+     */
+    void invokePriVarSwitch_(SolutionVector& uCurrentIter, std::true_type)
+    {
+        // update the variable switch (returns true if the pri vars at at least one dof were switched)
+        // for disabled grid variable caching
+        const auto& fvGridGeometry = assembler_->fvGridGeometry();
+        const auto& problem = assembler_->problem();
+        auto& gridVariables = assembler_->gridVariables();
+
+        // invoke the primary variable switch
+        priVarsSwitchedInLastIteration_ = priVarSwitch_->update(uCurrentIter, gridVariables,
+                                                                problem, fvGridGeometry);
+
+        if (priVarsSwitchedInLastIteration_)
+        {
+            for (const auto& element : elements(fvGridGeometry.gridView()))
+            {
+                // if the volume variables are cached globally, we need to update those where the primary variables have been switched
+                priVarSwitch_->updateSwitchedVolVars(problem, element, fvGridGeometry, gridVariables, uCurrentIter);
+
+                // if the flux variables are cached globally, we need to update those where the primary variables have been switched
+                priVarSwitch_->updateSwitchedFluxVarsCache(problem, element, fvGridGeometry, gridVariables, uCurrentIter);
+            }
+        }
+    }
+
     //! optimal number of iterations we want to achieve
     int targetSteps_;
     //! maximum number of iterations we do before giving up
@@ -786,7 +863,7 @@ private:
     //! assembleLinearSystem_ for assemblers that support partial reassembly
     template<class A>
     auto assembleLinearSystem_(const A& assembler, const SolutionVector& uCurrentIter)
-    -> typename std::enable_if_t<decltype(isValid(supportsPartialReassembly())(assembler))::value, void>
+    -> typename std::enable_if_t<decltype(isValid(Detail::supportsPartialReassembly())(assembler))::value, void>
     {
         assembler_->assembleJacobianAndResidual(uCurrentIter, partialReassembler_.get());
     }
@@ -794,7 +871,7 @@ private:
     //! assembleLinearSystem_ for assemblers that don't support partial reassembly
     template<class A>
     auto assembleLinearSystem_(const A& assembler, const SolutionVector& uCurrentIter)
-    -> typename std::enable_if_t<!decltype(isValid(supportsPartialReassembly())(assembler))::value, void>
+    -> typename std::enable_if_t<!decltype(isValid(Detail::supportsPartialReassembly())(assembler))::value, void>
     {
         assembler_->assembleJacobianAndResidual(uCurrentIter);
     }
@@ -1155,6 +1232,11 @@ private:
     std::size_t totalWastedIter_ = 0; //! Newton steps in solves that didn't converge
     std::size_t totalSucceededIter_ = 0; //! Newton steps in solves that converged
     std::size_t numConverged_ = 0; //! total number of converged solves
+
+    //! the class handling the primary variable switch
+    std::unique_ptr<PrimaryVariableSwitch> priVarSwitch_;
+    //! if we switched primary variables in the last iteration
+    bool priVarsSwitchedInLastIteration_ = false;
 };
 
 } // end namespace Dumux
diff --git a/dumux/nonlinear/privarswitchnewtonsolver.hh b/dumux/nonlinear/privarswitchnewtonsolver.hh
deleted file mode 100644
index dacf12d038f162442bb39bfac5d7ef0b5e0d0675..0000000000000000000000000000000000000000
--- a/dumux/nonlinear/privarswitchnewtonsolver.hh
+++ /dev/null
@@ -1,135 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/****************************************************************************
- *   See the file COPYING for full copying permissions.                      *
- *                                                                           *
- *   This program is free software: you can redistribute it and/or modify    *
- *   it under the terms of the GNU General Public License as published by    *
- *   the Free Software Foundation, either version 2 of the License, or       *
- *   (at your option) any later version.                                     *
- *                                                                           *
- *   This program is distributed in the hope that it will be useful,         *
- *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
- *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
- *   GNU General Public License for more details.                            *
- *                                                                           *
- *   You should have received a copy of the GNU General Public License       *
- *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
- *****************************************************************************/
-/*!
- * \file
- * \ingroup PorousmediumCompositional
- * \copydoc Dumux::PriVarSwitchNewtonSolver
- */
-#ifndef DUMUX_PRIVARSWITCH_NEWTON_SOLVER_HH
-#define DUMUX_PRIVARSWITCH_NEWTON_SOLVER_HH
-
-#include <memory>
-
-#include <dumux/common/parameters.hh>
-#include <dumux/discretization/method.hh>
-#include <dumux/discretization/elementsolution.hh>
-#include <dumux/nonlinear/newtonsolver.hh>
-
-namespace Dumux {
-
-/*!
- * \ingroup PorousmediumCompositional
- * \brief A newton solver that handles primary variable switches
- */
-template <class Assembler, class LinearSolver, class PrimaryVariableSwitch>
-class PriVarSwitchNewtonSolver : public NewtonSolver<Assembler, LinearSolver>
-{
-    using Scalar =  typename Assembler::Scalar;
-    using ParentType = NewtonSolver<Assembler, LinearSolver>;
-    using SolutionVector = typename Assembler::ResidualType;
-
-public:
-    using ParentType::ParentType;
-
-    /*!
-     * \brief Returns true if the error of the solution is below the
-     *        tolerance.
-     */
-    bool newtonConverged() const override
-    {
-        if (switchedInLastIteration_)
-            return false;
-
-        return ParentType::newtonConverged();
-    }
-
-    /*!
-     *
-     * \brief Called before the Newton method is applied to an
-     *        non-linear system of equations.
-     *
-     * \param u The initial solution
-     */
-    void newtonBegin(const SolutionVector &u) override
-    {
-        ParentType::newtonBegin(u);
-        const int verbosity = getParamFromGroup<int>(this->paramGroup(), "PrimaryVariableSwitch.Verbosity", 1);
-        priVarSwitch_ = std::make_unique<PrimaryVariableSwitch>(u.size(), verbosity);
-    }
-
-    /*!
-     * \brief Indicates that one Newton iteration was finished.
-     *
-     * \param assembler The jacobian assembler
-     * \param uCurrentIter The solution after the current Newton iteration
-     * \param uLastIter The solution at the beginning of the current Newton iteration
-     */
-    void newtonEndStep(SolutionVector &uCurrentIter,
-                       const SolutionVector &uLastIter) override
-    {
-        ParentType::newtonEndStep(uCurrentIter, uLastIter);
-
-        // update the variable switch (returns true if the pri vars at at least one dof were switched)
-        // for disabled grid variable caching
-        auto& assembler = this->assembler();
-        const auto& fvGridGeometry = assembler.fvGridGeometry();
-        const auto& problem = assembler.problem();
-        auto& gridVariables = assembler.gridVariables();
-
-        // invoke the primary variable switch
-        switchedInLastIteration_ = priVarSwitch_->update(uCurrentIter, gridVariables,
-                                                         problem, fvGridGeometry);
-
-        if(switchedInLastIteration_)
-        {
-            for (const auto& element : elements(fvGridGeometry.gridView()))
-            {
-                // if the volume variables are cached globally, we need to update those where the primary variables have been switched
-                priVarSwitch_->updateSwitchedVolVars(problem, element, fvGridGeometry, gridVariables, uCurrentIter);
-
-                // if the flux variables are cached globally, we need to update those where the primary variables have been switched
-                priVarSwitch_->updateSwitchedFluxVarsCache(problem, element, fvGridGeometry, gridVariables, uCurrentIter);
-            }
-        }
-    }
-
-    /*!
-     * \brief Called if the Newton method ended
-     *        (not known yet if we failed or succeeded)
-     */
-    void newtonEnd() override
-    {
-        ParentType::newtonEnd();
-
-        // in any way, we have to reset the switch flag
-        switchedInLastIteration_ = false;
-        // free some memory
-        priVarSwitch_.release();
-    }
-
-private:
-    //! the class handling the primary variable switch
-    std::unique_ptr<PrimaryVariableSwitch> priVarSwitch_;
-    //! if we switched primary variables in the last iteration
-    bool switchedInLastIteration_ = false;
-};
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/porousmediumflow/2p1c/model.hh b/dumux/porousmediumflow/2p1c/model.hh
index 4321cac320e7d06ae21186dad157545132a5cf5f..354f6f646292e7e90c0d0bfe78bfce9f59633155 100644
--- a/dumux/porousmediumflow/2p1c/model.hh
+++ b/dumux/porousmediumflow/2p1c/model.hh
@@ -76,7 +76,6 @@
 #include "localresidual.hh"
 #include "indices.hh"
 #include "volumevariables.hh"
-#include "primaryvariableswitch.hh"
 
 namespace Dumux {
 
@@ -191,10 +190,6 @@ public:
     using type = TwoPOneCVolumeVariables<Traits>;
 };
 
-//! The primary variable switch for the 2p1cni model.
-template<class TypeTag>
-struct PrimaryVariableSwitch<TypeTag, TTag::TwoPOneCNI> { using type = TwoPOneCPrimaryVariableSwitch; };
-
 //! The primary variables vector for the 2p1cni model.
 template<class TypeTag>
 struct PrimaryVariables<TypeTag, TTag::TwoPOneCNI>
diff --git a/dumux/porousmediumflow/2p1c/volumevariables.hh b/dumux/porousmediumflow/2p1c/volumevariables.hh
index 660029e45f0ef1b12d7b211f6d905eb68bb17064..230c9dccadaff363d25e74dfc2da44b877f0deaa 100644
--- a/dumux/porousmediumflow/2p1c/volumevariables.hh
+++ b/dumux/porousmediumflow/2p1c/volumevariables.hh
@@ -34,6 +34,8 @@
 #include <dumux/porousmediumflow/2p/formulation.hh>
 #include <dumux/material/solidstates/updatesolidvolumefractions.hh>
 
+#include "primaryvariableswitch.hh"
+
 namespace Dumux {
 
 /*!
@@ -91,6 +93,8 @@ public:
     using SolidState = typename Traits::SolidState;
     //! export type of solid system
     using SolidSystem = typename Traits::SolidSystem;
+    //! export the primary variable switch
+    using PrimaryVariableSwitch = TwoPOneCPrimaryVariableSwitch;
 
     //! return the two-phase formulation used here
     static constexpr TwoPFormulation priVarFormulation() { return formulation; }
diff --git a/dumux/porousmediumflow/2p2c/volumevariables.hh b/dumux/porousmediumflow/2p2c/volumevariables.hh
index 39ce7926ac8d190b46a5cbc89861a000c2ced5a0..21113f458037babafe9dadd6a1d317a628f01c30 100644
--- a/dumux/porousmediumflow/2p2c/volumevariables.hh
+++ b/dumux/porousmediumflow/2p2c/volumevariables.hh
@@ -34,6 +34,7 @@
 #include <dumux/porousmediumflow/nonisothermal/volumevariables.hh>
 #include <dumux/porousmediumflow/2p/formulation.hh>
 #include <dumux/material/solidstates/updatesolidvolumefractions.hh>
+#include <dumux/porousmediumflow/2pnc/primaryvariableswitch.hh>
 
 namespace Dumux {
 
@@ -93,6 +94,8 @@ public:
     using SolidState = typename Traits::SolidState;
     //! export type of solid system
     using SolidSystem = typename Traits::SolidSystem;
+    //! export the primary variable switch
+    using PrimaryVariableSwitch = TwoPNCPrimaryVariableSwitch;
 
     //! return whether moles or masses are balanced
     static constexpr bool useMoles() { return ModelTraits::useMoles(); }
diff --git a/dumux/porousmediumflow/2pnc/model.hh b/dumux/porousmediumflow/2pnc/model.hh
index f57efbf0d1ca2f164931db321d5a847cb3ca0c45..4dbd3baee8d540e38cc3fd9f20fd21804568451c 100644
--- a/dumux/porousmediumflow/2pnc/model.hh
+++ b/dumux/porousmediumflow/2pnc/model.hh
@@ -102,7 +102,6 @@
 #include <dumux/porousmediumflow/2p/formulation.hh>
 
 #include "volumevariables.hh"
-#include "primaryvariableswitch.hh"
 #include "iofields.hh"
 #include "indices.hh"
 
@@ -182,9 +181,6 @@ public:
     using type = SwitchablePrimaryVariables<PrimaryVariablesVector, int>;
 };
 
-template<class TypeTag>
-struct PrimaryVariableSwitch<TypeTag, TTag::TwoPNC> { using type = TwoPNCPrimaryVariableSwitch<TypeTag>; };         //!< The primary variable switch for the 2pnc model
-
 //! Set the volume variables property
 template<class TypeTag>
 struct VolumeVariables<TypeTag, TTag::TwoPNC>
diff --git a/dumux/porousmediumflow/2pnc/primaryvariableswitch.hh b/dumux/porousmediumflow/2pnc/primaryvariableswitch.hh
index 373af1a50c88dab1551151ce6f4a0198f6dd2e1d..3042f017076a4480b4c873ba0a9e4a33f6792b5b 100644
--- a/dumux/porousmediumflow/2pnc/primaryvariableswitch.hh
+++ b/dumux/porousmediumflow/2pnc/primaryvariableswitch.hh
@@ -35,11 +35,10 @@ namespace Dumux {
  * \ingroup TwoPNCModel
  * \brief The primary variable switch controlling the phase presence state variable
  */
-template<class TypeTag>
 class TwoPNCPrimaryVariableSwitch
-: public PrimaryVariableSwitch<TwoPNCPrimaryVariableSwitch<TypeTag>>
+: public PrimaryVariableSwitch<TwoPNCPrimaryVariableSwitch>
 {
-    using ParentType = PrimaryVariableSwitch<TwoPNCPrimaryVariableSwitch<TypeTag>>;
+    using ParentType = PrimaryVariableSwitch<TwoPNCPrimaryVariableSwitch>;
     friend ParentType;
 public:
     using ParentType::ParentType;
diff --git a/dumux/porousmediumflow/2pnc/volumevariables.hh b/dumux/porousmediumflow/2pnc/volumevariables.hh
index c0aa14e13b2035ef358c41dbc96a1f8e6aeca60f..25cebb9d811be45e63f1e90f1e8dc49b9d9858d0 100644
--- a/dumux/porousmediumflow/2pnc/volumevariables.hh
+++ b/dumux/porousmediumflow/2pnc/volumevariables.hh
@@ -41,6 +41,8 @@
 
 #include <dumux/porousmediumflow/2p/formulation.hh>
 
+#include "primaryvariableswitch.hh"
+
 namespace Dumux {
 
 /*!
@@ -98,6 +100,8 @@ public:
     using SolidState = typename Traits::SolidState;
     //! export type of solid system
     using SolidSystem = typename Traits::SolidSystem;
+    //! export the primary variable switch
+    using PrimaryVariableSwitch = TwoPNCPrimaryVariableSwitch;
 
     //! return whether moles or masses are balanced
     static constexpr bool useMoles() { return Traits::ModelTraits::useMoles(); }
diff --git a/dumux/porousmediumflow/3p3c/model.hh b/dumux/porousmediumflow/3p3c/model.hh
index 236e5da2bf9f02cab7bae4ca1762a69f7822e11c..cb80c0a7bf0815db70b79f0ec969fac13eb43ff7 100644
--- a/dumux/porousmediumflow/3p3c/model.hh
+++ b/dumux/porousmediumflow/3p3c/model.hh
@@ -96,7 +96,6 @@
 #include "indices.hh"
 #include "volumevariables.hh"
 #include "iofields.hh"
-#include "primaryvariableswitch.hh"
 #include "localresidual.hh"
 
 namespace Dumux {
@@ -198,10 +197,6 @@ struct FluidState<TypeTag, TTag::ThreePThreeC>{
 template<class TypeTag>
 struct LocalResidual<TypeTag, TTag::ThreePThreeC> { using type = ThreePThreeCLocalResidual<TypeTag>; };
 
-//! The primary variable switch for the 3p3c model
-template<class TypeTag>
-struct PrimaryVariableSwitch<TypeTag, TTag::ThreePThreeC> { using type = ThreePThreeCPrimaryVariableSwitch; };
-
 //! The primary variables vector for the 3p3c model
 template<class TypeTag>
 struct PrimaryVariables<TypeTag, TTag::ThreePThreeC>
diff --git a/dumux/porousmediumflow/3p3c/volumevariables.hh b/dumux/porousmediumflow/3p3c/volumevariables.hh
index 63099ebc933f318910eba5399117950f998b9281..190ea7b123b4f67fc305bb26294595850c061551 100644
--- a/dumux/porousmediumflow/3p3c/volumevariables.hh
+++ b/dumux/porousmediumflow/3p3c/volumevariables.hh
@@ -33,6 +33,8 @@
 #include <dumux/porousmediumflow/nonisothermal/volumevariables.hh>
 #include <dumux/material/solidstates/updatesolidvolumefractions.hh>
 
+#include "primaryvariableswitch.hh"
+
 namespace Dumux {
 
 /*!
@@ -93,6 +95,8 @@ public:
     using SolidState = typename Traits::SolidState;
     //! export type of solid system
     using SolidSystem = typename Traits::SolidSystem;
+    //! export the primary variable switch
+    using PrimaryVariableSwitch = ThreePThreeCPrimaryVariableSwitch;
 
     /*!
      * \brief Update all quantities for a given control volume
diff --git a/dumux/porousmediumflow/3pwateroil/model.hh b/dumux/porousmediumflow/3pwateroil/model.hh
index b9154dfd262bd086a41a332650e0bbd72243197a..bf5533536d507701489270a2f13670306be58813 100644
--- a/dumux/porousmediumflow/3pwateroil/model.hh
+++ b/dumux/porousmediumflow/3pwateroil/model.hh
@@ -88,7 +88,6 @@
 #include "model.hh"
 #include "volumevariables.hh"
 #include "localresidual.hh"
-#include "primaryvariableswitch.hh"
 #include "iofields.hh"
 
 namespace Dumux {
@@ -181,10 +180,6 @@ struct LocalResidual<TypeTag, TTag::ThreePWaterOilNI> { using type = ThreePWater
 template<class TypeTag>
 struct ReplaceCompEqIdx<TypeTag, TTag::ThreePWaterOilNI> { static constexpr int value = GetPropType<TypeTag, Properties::ModelTraits>::numComponents(); };
 
-//! The primary variable switch for the 3p3c model
-template<class TypeTag>
-struct PrimaryVariableSwitch<TypeTag, TTag::ThreePWaterOilNI> { using type = ThreePWaterOilPrimaryVariableSwitch<TypeTag>; };
-
 //! The primary variables vector for the 3p3c model
 template<class TypeTag>
 struct PrimaryVariables<TypeTag, TTag::ThreePWaterOilNI>
diff --git a/dumux/porousmediumflow/3pwateroil/primaryvariableswitch.hh b/dumux/porousmediumflow/3pwateroil/primaryvariableswitch.hh
index de858bd0c5a1b50461be1f784f38b802116fb335..e3c889722456229d6a2c0414bfc22709883a8d92 100644
--- a/dumux/porousmediumflow/3pwateroil/primaryvariableswitch.hh
+++ b/dumux/porousmediumflow/3pwateroil/primaryvariableswitch.hh
@@ -33,92 +33,65 @@ namespace Dumux {
  * \ingroup ThreePWaterOilModel
  * \brief The primary variable switch controlling the phase presence state variable
  */
-template<class TypeTag>
 class ThreePWaterOilPrimaryVariableSwitch
-: public PrimaryVariableSwitch<ThreePWaterOilPrimaryVariableSwitch<TypeTag>>
+: public PrimaryVariableSwitch<ThreePWaterOilPrimaryVariableSwitch>
 {
-    using ParentType = PrimaryVariableSwitch<ThreePWaterOilPrimaryVariableSwitch<TypeTag>>;
+    using ParentType = PrimaryVariableSwitch<ThreePWaterOilPrimaryVariableSwitch>;
     friend ParentType;
 
-    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
-    using GridView = GetPropType<TypeTag, Properties::GridView>;
-    using IndexType = typename GridView::IndexSet::IndexType;
-    using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimensionworld>;
-
-    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
-    using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
-    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
-    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
-
-    enum {
-        switch1Idx = Indices::switch1Idx,
-        switch2Idx = Indices::switch2Idx,
-        pressureIdx = Indices::pressureIdx,
-
-        wPhaseIdx = FluidSystem::wPhaseIdx,
-        nPhaseIdx = FluidSystem::nPhaseIdx,
-        gPhaseIdx = FluidSystem::gPhaseIdx,
-
-        wCompIdx = FluidSystem::wCompIdx,
-        nCompIdx = FluidSystem::nCompIdx,
-
-        threePhases = Indices::threePhases,
-        wPhaseOnly  = Indices::wPhaseOnly,
-        gnPhaseOnly = Indices::gnPhaseOnly,
-        wnPhaseOnly = Indices::wnPhaseOnly,
-        gPhaseOnly  = Indices::gPhaseOnly,
-        wgPhaseOnly = Indices::wgPhaseOnly
-    };
-
 public:
     using ParentType::ParentType;
 
 protected:
 
     // perform variable switch at a degree of freedom location
-    bool update_(PrimaryVariables& priVars,
+    template<class VolumeVariables, class GlobalPosition>
+    bool update_(typename VolumeVariables::PrimaryVariables& priVars,
                  const VolumeVariables& volVars,
-                 IndexType dofIdxGlobal,
+                 std::size_t dofIdxGlobal,
                  const GlobalPosition& globalPos)
     {
+        using PrimaryVariables = typename VolumeVariables::PrimaryVariables;
+        using Scalar = typename PrimaryVariables::value_type;
+        using Indices = typename VolumeVariables::Indices;
+        using FluidSystem = typename VolumeVariables::FluidSystem;
+
         // evaluate primary variable switch
         bool wouldSwitch = false;
         auto phasePresence = priVars.state();
         int newPhasePresence = phasePresence;
 
-        bool onlyGasPhaseCanDisappear = getPropValue<TypeTag, Properties::OnlyGasPhaseCanDisappear>();
-
-        if(onlyGasPhaseCanDisappear)
+        if(VolumeVariables::onlyGasPhaseCanDisappear())
         {
             // check if a primary var switch is necessary
-            if (phasePresence == threePhases)
+            if (phasePresence == Indices::threePhases)
             {
                 Scalar Smin = 0;
                 if (this->wasSwitched_[dofIdxGlobal])
                     Smin = -0.01;
 
-                if (volVars.saturation(gPhaseIdx) <= Smin)
+                if (volVars.saturation(FluidSystem::gPhaseIdx) <= Smin)
                 {
                     wouldSwitch = true;
                     // gas phase disappears
                     if (this->verbosity() > 1)
                         std::cout << "Gas phase disappears at dof " << dofIdxGlobal
                                 << ", coordinates: " << globalPos << ", sg: "
-                                << volVars.saturation(gPhaseIdx) << std::endl;
-                    newPhasePresence = wnPhaseOnly;
+                                << volVars.saturation(FluidSystem::gPhaseIdx) << std::endl;
+                    newPhasePresence = FluidSystem::gPhaseIdx;
 
-                    priVars[pressureIdx] = volVars.fluidState().pressure(wPhaseIdx);
-                    priVars[switch1Idx] = volVars.fluidState().temperature();
+                    priVars[Indices::pressureIdx] = volVars.fluidState().pressure(FluidSystem::wPhaseIdx);
+                    priVars[Indices::switch1Idx] = volVars.fluidState().temperature();
                 }
             }
-            else if (phasePresence == wnPhaseOnly)
+            else if (phasePresence == FluidSystem::gPhaseIdx)
             {
                 bool gasFlag = 0;
 
                 // calculate fractions of the partial pressures in the
                 // hypothetical gas phase
-                Scalar xwg = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
-                Scalar xng = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
+                Scalar xwg = volVars.fluidState().moleFraction(FluidSystem::gPhaseIdx, FluidSystem::wCompIdx);
+                Scalar xng = volVars.fluidState().moleFraction(FluidSystem::gPhaseIdx, FluidSystem::nCompIdx);
 
                 Scalar xgMax = 1.0;
                 if (xwg + xng > xgMax)
@@ -140,9 +113,9 @@ protected:
 
                 if (gasFlag == 1)
                 {
-                    newPhasePresence = threePhases;
-                    priVars[pressureIdx] = volVars.pressure(gPhaseIdx);
-                    priVars[switch1Idx] = volVars.saturation(wPhaseIdx);
+                    newPhasePresence = Indices::threePhases;
+                    priVars[Indices::pressureIdx] = volVars.pressure(FluidSystem::gPhaseIdx);
+                    priVars[Indices::switch1Idx] = volVars.saturation(FluidSystem::wPhaseIdx);
                 }
             }
         }
@@ -150,59 +123,59 @@ protected:
         else
         {
             // check if a primary var switch is necessary
-            if (phasePresence == threePhases)
+            if (phasePresence == Indices::threePhases)
             {
                 Scalar Smin = 0;
                 if (this->wasSwitched_[dofIdxGlobal])
                     Smin = -0.01;
 
-                if (volVars.saturation(gPhaseIdx) <= Smin)
+                if (volVars.saturation(FluidSystem::gPhaseIdx) <= Smin)
                 {
                     wouldSwitch = true;
                     // gas phase disappears
                     if (this->verbosity() > 1)
                         std::cout << "Gas phase disappears at dof " << dofIdxGlobal
                                 << ", coordinates: " << globalPos << ", sg: "
-                                << volVars.saturation(gPhaseIdx) << std::endl;
-                    newPhasePresence = wnPhaseOnly;
+                                << volVars.saturation(FluidSystem::gPhaseIdx) << std::endl;
+                    newPhasePresence = FluidSystem::gPhaseIdx;
 
-                    priVars[pressureIdx] = volVars.fluidState().pressure(wPhaseIdx);
-                    priVars[switch1Idx] = volVars.fluidState().temperature();
+                    priVars[Indices::pressureIdx] = volVars.fluidState().pressure(FluidSystem::wPhaseIdx);
+                    priVars[Indices::switch1Idx] = volVars.fluidState().temperature();
                 }
-                else if (volVars.saturation(wPhaseIdx) <= Smin)
+                else if (volVars.saturation(FluidSystem::wPhaseIdx) <= Smin)
                 {
                     wouldSwitch = true;
                     // water phase disappears
                     if (this->verbosity() > 1)
                         std::cout << "Water phase disappears at dof " << dofIdxGlobal
                                 << ", coordinates: " << globalPos << ", sw: "
-                                << volVars.saturation(wPhaseIdx) << std::endl;
-                    newPhasePresence = gnPhaseOnly;
+                                << volVars.saturation(FluidSystem::wPhaseIdx) << std::endl;
+                    newPhasePresence = Indices::gnPhaseOnly;
 
-                    priVars[switch1Idx] = volVars.fluidState().saturation(nPhaseIdx);
-                    priVars[switch2Idx] = volVars.fluidState().moleFraction(nPhaseIdx, wCompIdx);
+                    priVars[Indices::switch1Idx] = volVars.fluidState().saturation(FluidSystem::nPhaseIdx);
+                    priVars[Indices::switch2Idx] = volVars.fluidState().moleFraction(FluidSystem::nPhaseIdx, FluidSystem::wCompIdx);
                 }
-                else if (volVars.saturation(nPhaseIdx) <= Smin)
+                else if (volVars.saturation(FluidSystem::nPhaseIdx) <= Smin)
                 {
                     wouldSwitch = true;
                     // NAPL phase disappears
                     if (this->verbosity() > 1)
                         std::cout << "NAPL phase disappears at dof " << dofIdxGlobal
                                 << ", coordinates: " << globalPos << ", sn: "
-                                << volVars.saturation(nPhaseIdx) << std::endl;
-                    newPhasePresence = wgPhaseOnly;
+                                << volVars.saturation(FluidSystem::nPhaseIdx) << std::endl;
+                    newPhasePresence = Indices::wgPhaseOnly;
 
-                    priVars[switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
+                    priVars[Indices::switch2Idx] = volVars.fluidState().moleFraction(FluidSystem::wPhaseIdx, FluidSystem::nCompIdx);
                 }
             }
-            else if (phasePresence == wPhaseOnly)
+            else if (phasePresence == Indices::wPhaseOnly)
             {
                 bool gasFlag = 0;
                 bool nonwettingFlag = 0;
                 // calculate fractions of the partial pressures in the
                 // hypothetical gas phase
-                Scalar xwg = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
-                Scalar xng = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
+                Scalar xwg = volVars.fluidState().moleFraction(FluidSystem::gPhaseIdx, FluidSystem::wCompIdx);
+                Scalar xng = volVars.fluidState().moleFraction(FluidSystem::gPhaseIdx, FluidSystem::nCompIdx);
 
                 Scalar xgMax = 1.0;
                 if (xwg + xng > xgMax)
@@ -223,7 +196,7 @@ protected:
                 }
 
                 // calculate fractions in the hypothetical NAPL phase
-                Scalar xnn = volVars.fluidState().moleFraction(nPhaseIdx, nCompIdx);
+                Scalar xnn = volVars.fluidState().moleFraction(FluidSystem::nPhaseIdx, FluidSystem::nCompIdx);
                 /* take care:
                 for xnn in case wPhaseOnly we compute xnn=henry_mesitylene*x1w,
                 where a hypothetical gas pressure is assumed for the Henry
@@ -251,23 +224,23 @@ protected:
 
                 if ((gasFlag == 1) && (nonwettingFlag == 0))
                 {
-                    newPhasePresence = wgPhaseOnly;
-                    priVars[switch1Idx] = 0.9999;
-                    priVars[pressureIdx] = volVars.fluidState().pressure(gPhaseIdx);
+                    newPhasePresence = Indices::wgPhaseOnly;
+                    priVars[Indices::switch1Idx] = 0.9999;
+                    priVars[Indices::pressureIdx] = volVars.fluidState().pressure(FluidSystem::gPhaseIdx);
                 }
                 else if ((gasFlag == 1) && (nonwettingFlag == 1))
                 {
-                    newPhasePresence = threePhases;
-                    priVars[pressureIdx] = volVars.fluidState().pressure(gPhaseIdx);
-                    priVars[switch1Idx] = 0.999;
+                    newPhasePresence = Indices::threePhases;
+                    priVars[Indices::pressureIdx] = volVars.fluidState().pressure(FluidSystem::gPhaseIdx);
+                    priVars[Indices::switch1Idx] = 0.999;
                 }
                 else if ((gasFlag == 0) && (nonwettingFlag == 1))
                 {
-                    newPhasePresence = wnPhaseOnly;
-                    priVars[switch2Idx] = 0.0001;
+                    newPhasePresence = Indices::wnPhaseOnly;
+                    priVars[Indices::switch2Idx] = 0.0001;
                 }
             }
-            else if (phasePresence == gnPhaseOnly)
+            else if (phasePresence == Indices::gnPhaseOnly)
             {
                 bool nonwettingFlag = 0;
                 bool wettingFlag = 0;
@@ -276,20 +249,20 @@ protected:
                 if (this->wasSwitched_[dofIdxGlobal])
                     Smin = -0.01;
 
-                if (volVars.saturation(nPhaseIdx) <= Smin)
+                if (volVars.saturation(FluidSystem::nPhaseIdx) <= Smin)
                 {
                     wouldSwitch = true;
                     // NAPL phase disappears
                     if (this->verbosity() > 1)
                         std::cout << "NAPL phase disappears at dof " << dofIdxGlobal
                                 << ", coordinates: " << globalPos << ", sn: "
-                                << volVars.saturation(nPhaseIdx) << std::endl;
+                                << volVars.saturation(FluidSystem::nPhaseIdx) << std::endl;
                     nonwettingFlag = 1;
                 }
 
 
                 // calculate fractions of the hypothetical water phase
-                Scalar xww = volVars.fluidState().moleFraction(wPhaseIdx, wCompIdx);
+                Scalar xww = volVars.fluidState().moleFraction(FluidSystem::wPhaseIdx, FluidSystem::wCompIdx);
                 /*
                 take care:, xww, if no water is present, then take xww=xwg*pg/pwsat .
                 If this is larger than 1, then water appears
@@ -314,24 +287,24 @@ protected:
 
                 if ((wettingFlag == 1) && (nonwettingFlag == 0))
                 {
-                    newPhasePresence = threePhases;
-                    priVars[switch1Idx] = 0.0001;
-                    priVars[switch2Idx] = volVars.saturation(nPhaseIdx);
+                    newPhasePresence = Indices::threePhases;
+                    priVars[Indices::switch1Idx] = 0.0001;
+                    priVars[Indices::switch2Idx] = volVars.saturation(FluidSystem::nPhaseIdx);
                 }
                 else if ((wettingFlag == 1) && (nonwettingFlag == 1))
                 {
-                    newPhasePresence = wgPhaseOnly;
-                    priVars[switch1Idx] = 0.0001;
-                    priVars[switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
+                    newPhasePresence = Indices::wgPhaseOnly;
+                    priVars[Indices::switch1Idx] = 0.0001;
+                    priVars[Indices::switch2Idx] = volVars.fluidState().moleFraction(FluidSystem::wPhaseIdx, FluidSystem::nCompIdx);
                 }
                 else if ((wettingFlag == 0) && (nonwettingFlag == 1))
                 {
-                    newPhasePresence = gPhaseOnly;
-                    priVars[switch1Idx] = volVars.fluidState().temperature();
-                    priVars[switch2Idx] = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
+                    newPhasePresence = Indices::gPhaseOnly;
+                    priVars[Indices::switch1Idx] = volVars.fluidState().temperature();
+                    priVars[Indices::switch2Idx] = volVars.fluidState().moleFraction(FluidSystem::gPhaseIdx, FluidSystem::nCompIdx);
                 }
             }
-            else if (phasePresence == wnPhaseOnly)
+            else if (phasePresence == Indices::wnPhaseOnly)
             {
                 bool nonwettingFlag = 0;
                 bool gasFlag = 0;
@@ -340,21 +313,21 @@ protected:
                 if (this->wasSwitched_[dofIdxGlobal])
                     Smin = -0.01;
 
-                if (volVars.saturation(nPhaseIdx) <= Smin)
+                if (volVars.saturation(FluidSystem::nPhaseIdx) <= Smin)
                 {
                     wouldSwitch = true;
                     // NAPL phase disappears
                     if (this->verbosity() > 1)
                         std::cout << "NAPL phase disappears at dof " << dofIdxGlobal
                                 << ", coordinates: " << globalPos << ", sn: "
-                                << volVars.saturation(nPhaseIdx) << std::endl;
+                                << volVars.saturation(FluidSystem::nPhaseIdx) << std::endl;
                     nonwettingFlag = 1;
                 }
 
                 // calculate fractions of the partial pressures in the
                 // hypothetical gas phase
-                Scalar xwg = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
-                Scalar xng = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
+                Scalar xwg = volVars.fluidState().moleFraction(FluidSystem::gPhaseIdx, FluidSystem::wCompIdx);
+                Scalar xng = volVars.fluidState().moleFraction(FluidSystem::gPhaseIdx, FluidSystem::nCompIdx);
 
                 Scalar xgMax = 1.0;
                 if (xwg + xng > xgMax)
@@ -376,30 +349,30 @@ protected:
 
                 if ((gasFlag == 1) && (nonwettingFlag == 0))
                 {
-                    newPhasePresence = threePhases;
-                    priVars[pressureIdx] = volVars.pressure(gPhaseIdx);
-                    priVars[switch1Idx] = volVars.saturation(wPhaseIdx);
+                    newPhasePresence = Indices::threePhases;
+                    priVars[Indices::pressureIdx] = volVars.pressure(FluidSystem::gPhaseIdx);
+                    priVars[Indices::switch1Idx] = volVars.saturation(FluidSystem::wPhaseIdx);
                 }
                 else if ((gasFlag == 1) && (nonwettingFlag == 1))
                 {
-                    newPhasePresence = wgPhaseOnly;
-                    priVars[pressureIdx] = volVars.pressure(gPhaseIdx);
-                    priVars[switch1Idx] = volVars.temperature();
-                    priVars[switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
+                    newPhasePresence = Indices::wgPhaseOnly;
+                    priVars[Indices::pressureIdx] = volVars.pressure(FluidSystem::gPhaseIdx);
+                    priVars[Indices::switch1Idx] = volVars.temperature();
+                    priVars[Indices::switch2Idx] = volVars.fluidState().moleFraction(FluidSystem::wPhaseIdx, FluidSystem::nCompIdx);
                 }
                 else if ((gasFlag == 0) && (nonwettingFlag == 1))
                 {
-                    newPhasePresence = wPhaseOnly;
-                    priVars[switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
+                    newPhasePresence = Indices::wPhaseOnly;
+                    priVars[Indices::switch2Idx] = volVars.fluidState().moleFraction(FluidSystem::wPhaseIdx, FluidSystem::nCompIdx);
                 }
             }
-            else if (phasePresence == gPhaseOnly)
+            else if (phasePresence == Indices::gPhaseOnly)
             {
                 bool nonwettingFlag = 0;
                 bool wettingFlag = 0;
 
                 // calculate fractions in the hypothetical NAPL phase
-                Scalar xnn = volVars.fluidState().moleFraction(nPhaseIdx, nCompIdx);
+                Scalar xnn = volVars.fluidState().moleFraction(FluidSystem::nPhaseIdx, FluidSystem::nCompIdx);
                 /*
                 take care:, xnn, if no NAPL phase is there, take xnn=xng*pg/pcsat
                 if this is larger than 1, then NAPL appears
@@ -423,7 +396,7 @@ protected:
                     nonwettingFlag = 1;
                 }
                 // calculate fractions of the hypothetical water phase
-                Scalar xww = volVars.fluidState().moleFraction(wPhaseIdx, wCompIdx);
+                Scalar xww = volVars.fluidState().moleFraction(FluidSystem::wPhaseIdx, FluidSystem::wCompIdx);
                 /*
                 take care:, xww, if no water is present, then take xww=xwg*pg/pwsat .
                 If this is larger than 1, then water appears
@@ -447,31 +420,31 @@ protected:
                 }
                 if ((wettingFlag == 1) && (nonwettingFlag == 0))
                 {
-                    newPhasePresence = wgPhaseOnly;
-                    priVars[switch1Idx] = 0.0001;
-                    priVars[switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
+                    newPhasePresence = Indices::wgPhaseOnly;
+                    priVars[Indices::switch1Idx] = 0.0001;
+                    priVars[Indices::switch2Idx] = volVars.fluidState().moleFraction(FluidSystem::wPhaseIdx, FluidSystem::nCompIdx);
                 }
                 else if ((wettingFlag == 1) && (nonwettingFlag == 1))
                 {
-                    newPhasePresence = threePhases;
-                    priVars[switch1Idx] = 0.0001;
-                    priVars[switch2Idx] = 0.0001;
+                    newPhasePresence = Indices::threePhases;
+                    priVars[Indices::switch1Idx] = 0.0001;
+                    priVars[Indices::switch2Idx] = 0.0001;
                 }
                 else if ((wettingFlag == 0) && (nonwettingFlag == 1))
                 {
-                    newPhasePresence = gnPhaseOnly;
-                    priVars[switch2Idx]
-                        = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
+                    newPhasePresence = Indices::gnPhaseOnly;
+                    priVars[Indices::switch2Idx]
+                        = volVars.fluidState().moleFraction(FluidSystem::wPhaseIdx, FluidSystem::nCompIdx);
                 }
             }
-            else if (phasePresence == wgPhaseOnly)
+            else if (phasePresence == Indices::wgPhaseOnly)
             {
                 bool nonwettingFlag = 0;
                 bool gasFlag = 0;
                 bool wettingFlag = 0;
 
                 // get the fractions in the hypothetical NAPL phase
-                Scalar xnn = volVars.fluidState().moleFraction(nPhaseIdx, nCompIdx);
+                Scalar xnn = volVars.fluidState().moleFraction(FluidSystem::nPhaseIdx, FluidSystem::nCompIdx);
 
                 // take care: if the NAPL phase is not present, take
                 // xnn=xng*pg/pcsat if this is larger than 1, then NAPL
@@ -498,14 +471,14 @@ protected:
                 if (this->wasSwitched_[dofIdxGlobal])
                     Smin = -0.01;
 
-                if (volVars.saturation(gPhaseIdx) <= Smin)
+                if (volVars.saturation(FluidSystem::gPhaseIdx) <= Smin)
                 {
                     wouldSwitch = true;
                     // gas phase disappears
                     if (this->verbosity() > 1)
                         std::cout << "Gas phase disappears at dof " << dofIdxGlobal
                                 << ", coordinates: " << globalPos << ", sg: "
-                                << volVars.saturation(gPhaseIdx) << std::endl;
+                                << volVars.saturation(FluidSystem::gPhaseIdx) << std::endl;
                     gasFlag = 1;
                 }
 
@@ -513,43 +486,43 @@ protected:
                 if (this->wasSwitched_[dofIdxGlobal])
                     Smin = -0.01;
 
-                if (volVars.saturation(wPhaseIdx) <= Smin)
+                if (volVars.saturation(FluidSystem::wPhaseIdx) <= Smin)
                 {
                     wouldSwitch = true;
                     // gas phase disappears
                     if (this->verbosity() > 1)
                         std::cout << "Water phase disappears at dof " << dofIdxGlobal
                                 << ", coordinates: " << globalPos << ", sw: "
-                                << volVars.saturation(wPhaseIdx) << std::endl;
+                                << volVars.saturation(FluidSystem::wPhaseIdx) << std::endl;
                     wettingFlag = 1;
                 }
 
                 if ((gasFlag == 0) && (nonwettingFlag == 1) && (wettingFlag == 1))
                 {
-                    newPhasePresence = gnPhaseOnly;
-                    priVars[switch1Idx] = 0.0001;
-                    priVars[switch2Idx] = volVars.fluidState().moleFraction(nPhaseIdx, wCompIdx);
+                    newPhasePresence = Indices::gnPhaseOnly;
+                    priVars[Indices::switch1Idx] = 0.0001;
+                    priVars[Indices::switch2Idx] = volVars.fluidState().moleFraction(FluidSystem::nPhaseIdx, FluidSystem::wCompIdx);
     ;
                 }
                 else if ((gasFlag == 0) && (nonwettingFlag == 1) && (wettingFlag == 0))
                 {
-                    newPhasePresence = threePhases;
-                    priVars[switch2Idx] = 0.0001;
+                    newPhasePresence = Indices::threePhases;
+                    priVars[Indices::switch2Idx] = 0.0001;
                 }
                 else if ((gasFlag == 1) && (nonwettingFlag == 0) && (wettingFlag == 0))
                 {
-                    newPhasePresence = wPhaseOnly;
-                    priVars[pressureIdx] = volVars.fluidState().pressure(wPhaseIdx);
-                    priVars[switch1Idx] = volVars.fluidState().temperature();
-                    priVars[switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
+                    newPhasePresence = Indices::wPhaseOnly;
+                    priVars[Indices::pressureIdx] = volVars.fluidState().pressure(FluidSystem::wPhaseIdx);
+                    priVars[Indices::switch1Idx] = volVars.fluidState().temperature();
+                    priVars[Indices::switch2Idx] = volVars.fluidState().moleFraction(FluidSystem::wPhaseIdx, FluidSystem::nCompIdx);
                 }
                 else if ((gasFlag == 0) && (nonwettingFlag == 0) && (wettingFlag == 1))
                 {
-                    newPhasePresence = gPhaseOnly;
-                    priVars[switch1Idx]
+                    newPhasePresence = Indices::gPhaseOnly;
+                    priVars[Indices::switch1Idx]
                         = volVars.fluidState().temperature();
-                    priVars[switch2Idx]
-                        = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
+                    priVars[Indices::switch2Idx]
+                        = volVars.fluidState().moleFraction(FluidSystem::gPhaseIdx, FluidSystem::nCompIdx);
                 }
             }
         }
diff --git a/dumux/porousmediumflow/3pwateroil/volumevariables.hh b/dumux/porousmediumflow/3pwateroil/volumevariables.hh
index eae4cf8e82d604443098a47a8a720489dece1ce3..e30b252f28efb714f37c3177c598010808e9dec7 100644
--- a/dumux/porousmediumflow/3pwateroil/volumevariables.hh
+++ b/dumux/porousmediumflow/3pwateroil/volumevariables.hh
@@ -39,6 +39,8 @@
 #include <dumux/common/valgrind.hh>
 #include <dumux/common/exceptions.hh>
 
+#include "primaryvariableswitch.hh"
+
 namespace Dumux {
 
 /*!
@@ -55,7 +57,6 @@ class ThreePWaterOilVolumeVariables
     using EnergyVolVars = EnergyVolumeVariables<Traits, ThreePWaterOilVolumeVariables<Traits> >;
     using Scalar = typename Traits::PrimaryVariables::value_type;
     using ModelTraits = typename Traits::ModelTraits;
-    using Indices = typename ModelTraits::Indices;
     using FS = typename Traits::FluidSystem;
     static constexpr int numFluidComps = ParentType::numComponents();
 
@@ -69,19 +70,19 @@ class ThreePWaterOilVolumeVariables
         gPhaseIdx = FS::gPhaseIdx,
         nPhaseIdx = FS::nPhaseIdx,
 
-        switch1Idx = Indices::switch1Idx,
-        switch2Idx = Indices::switch2Idx,
-        pressureIdx = Indices::pressureIdx
+        switch1Idx = ModelTraits::Indices::switch1Idx,
+        switch2Idx = ModelTraits::Indices::switch2Idx,
+        pressureIdx = ModelTraits::Indices::pressureIdx
     };
 
     // present phases
     enum {
-        threePhases = Indices::threePhases,
-        wPhaseOnly  = Indices::wPhaseOnly,
-        gnPhaseOnly = Indices::gnPhaseOnly,
-        wnPhaseOnly = Indices::wnPhaseOnly,
-        gPhaseOnly  = Indices::gPhaseOnly,
-        wgPhaseOnly = Indices::wgPhaseOnly
+        threePhases = ModelTraits::Indices::threePhases,
+        wPhaseOnly  = ModelTraits::Indices::wPhaseOnly,
+        gnPhaseOnly = ModelTraits::Indices::gnPhaseOnly,
+        wnPhaseOnly = ModelTraits::Indices::wnPhaseOnly,
+        gPhaseOnly  = ModelTraits::Indices::gPhaseOnly,
+        wgPhaseOnly = ModelTraits::Indices::wgPhaseOnly
     };
 
 public:
@@ -89,10 +90,17 @@ public:
     using FluidState = typename Traits::FluidState;
     //! The type of the fluid system
     using FluidSystem = typename Traits::FluidSystem;
+    //! export the indices
+    using Indices = typename ModelTraits::Indices;
     //! export type of solid state
     using SolidState = typename Traits::SolidState;
     //! export type of solid system
     using SolidSystem = typename Traits::SolidSystem;
+    //! export the primary variable switch
+    using PrimaryVariableSwitch = ThreePWaterOilPrimaryVariableSwitch;
+    //! state if only the gas phase is allowed to disappear
+    static constexpr bool onlyGasPhaseCanDisappear()
+    { return Traits::ModelTraits::onlyGasPhaseCanDisappear(); }
 
     /*!
      * \copydoc ImplicitVolumeVariables::update
@@ -107,13 +115,11 @@ public:
         const auto& priVars = elemSol[scv.localDofIndex()];
         const auto phasePresence = priVars.state();
 
-        bool onlyGasPhaseCanDisappear = Traits::ModelTraits::onlyGasPhaseCanDisappear();
-
         // capillary pressure parameters
         using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
         const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
 
-        if(!onlyGasPhaseCanDisappear)
+        if(!onlyGasPhaseCanDisappear())
         {
             /* first the saturations */
             if (phasePresence == threePhases)
diff --git a/dumux/porousmediumflow/co2/model.hh b/dumux/porousmediumflow/co2/model.hh
index cb9f45dbc0a4db8cd9a7e38ac9eae59801170573..a57a6a0e3d4dd10e1296284b5f4590e47c66c14d 100644
--- a/dumux/porousmediumflow/co2/model.hh
+++ b/dumux/porousmediumflow/co2/model.hh
@@ -26,7 +26,6 @@
 
 #include <dumux/common/properties.hh>
 #include <dumux/porousmediumflow/2p2c/model.hh>
-#include "primaryvariableswitch.hh"
 #include "volumevariables.hh"
 
 /*!
@@ -80,12 +79,6 @@ struct TwoPTwoCCO2 { using InheritsFrom = std::tuple<TwoPTwoC>; };
 struct TwoPTwoCCO2NI { using InheritsFrom = std::tuple<TwoPTwoCNI>; };
 } // end namespace TTag
 
-//! the CO2 privarswitch and VolumeVariables properties
-template<class TypeTag>
-struct PrimaryVariableSwitch<TypeTag, TTag::TwoPTwoCCO2> { using type = TwoPTwoCCO2PrimaryVariableSwitch; };
-template<class TypeTag>
-struct PrimaryVariableSwitch<TypeTag, TTag::TwoPTwoCCO2NI> { using type = TwoPTwoCCO2PrimaryVariableSwitch; };
-
 //! the co2 volume variables use the same traits as the 2p2c model
 template<class TypeTag>
 struct VolumeVariables<TypeTag, TTag::TwoPTwoCCO2>
diff --git a/dumux/porousmediumflow/co2/volumevariables.hh b/dumux/porousmediumflow/co2/volumevariables.hh
index c17aab3f86179e992c94ed48fd6f14f1ec4bc8e5..1691be48d362be97a6a5c31e12b6a7b25cd759da 100644
--- a/dumux/porousmediumflow/co2/volumevariables.hh
+++ b/dumux/porousmediumflow/co2/volumevariables.hh
@@ -34,6 +34,8 @@
 #include <dumux/porousmediumflow/nonisothermal/volumevariables.hh>
 #include <dumux/material/solidstates/updatesolidvolumefractions.hh>
 
+#include "primaryvariableswitch.hh"
+
 namespace Dumux {
 
 /*!
@@ -91,6 +93,8 @@ public:
     using SolidState = typename Traits::SolidState;
     //! export type of solid system
     using SolidSystem = typename Traits::SolidSystem;
+    //! export the type of the primary variable switch
+    using PrimaryVariableSwitch = TwoPTwoCCO2PrimaryVariableSwitch;
 
 
     //! return whether moles or masses are balanced
diff --git a/dumux/porousmediumflow/compositional/primaryvariableswitch.hh b/dumux/porousmediumflow/compositional/primaryvariableswitch.hh
index daf8668dd4d706702b57576d087bef3b7aae4e1a..bfb99fc736c6de1fd04215f775364a9e6a1a759e 100644
--- a/dumux/porousmediumflow/compositional/primaryvariableswitch.hh
+++ b/dumux/porousmediumflow/compositional/primaryvariableswitch.hh
@@ -43,12 +43,11 @@ public:
     template<typename... Args>
     NoPrimaryVariableSwitch(Args&&...) {}
 
-    template<typename... Args> void init(Args&&...) {}
+    template<typename... Args> void reset(Args&&...) {}
     template<typename... Args> bool wasSwitched(Args&&...) const { return false; }
     template<typename... Args> bool update(Args&&...) { return false; }
     template<typename... Args> void updateSwitchedVolVars(Args&&...) {}
     template<typename... Args> void updateSwitchedFluxVarsCache(Args&&...) {}
-    template<typename... Args> bool update_(Args&&...) {return false; }
 };
 
 /*!
@@ -59,11 +58,9 @@ template<class Implementation>
 class PrimaryVariableSwitch
 {
 public:
-    PrimaryVariableSwitch(const std::size_t& numDofs, int verbosity = 1)
+    PrimaryVariableSwitch(int verbosity = 1)
     : verbosity_(verbosity)
-    {
-        wasSwitched_.resize(numDofs, false);
-    }
+    {}
 
     //! If the primary variables were recently switched
     bool wasSwitched(std::size_t dofIdxGlobal) const
@@ -71,6 +68,12 @@ public:
         return wasSwitched_[dofIdxGlobal];
     }
 
+    //! Reset all flags
+    void reset(const std::size_t numDofs)
+    {
+        wasSwitched_.resize(numDofs, false);
+    }
+
     /*!
      * \brief Update the variable switch / phase presence
      *
diff --git a/dumux/porousmediumflow/properties.hh b/dumux/porousmediumflow/properties.hh
index 54273367d9a5f4ddc6971418354a72d1bccee1db..30dbff42bf3c75ec28e424d43af0ef2cc8f3e36a 100644
--- a/dumux/porousmediumflow/properties.hh
+++ b/dumux/porousmediumflow/properties.hh
@@ -33,7 +33,6 @@
 #include <dumux/porousmediumflow/fluxvariables.hh>
 #include <dumux/porousmediumflow/fluxvariablescache.hh>
 #include <dumux/porousmediumflow/nonisothermal/localresidual.hh>
-#include <dumux/porousmediumflow/compositional/primaryvariableswitch.hh>
 #include <dumux/porousmediumflow/velocityoutput.hh>
 
 #include <dumux/flux/darcyslaw.hh>
@@ -93,10 +92,6 @@ struct VelocityOutput<TypeTag, TTag::PorousMediumFlow>
                                                 GetPropType<TypeTag, Properties::FluxVariables>>;
 };
 
-//! By default, we set an empty primary variables switch
-template<class TypeTag>
-struct PrimaryVariableSwitch<TypeTag, TTag::PorousMediumFlow> { using type = NoPrimaryVariableSwitch; };
-
 template<class TypeTag>
 struct EnableThermalNonEquilibrium<TypeTag, TTag::PorousMediumFlow> { static constexpr bool value = false; };
 
diff --git a/dumux/porousmediumflow/richards/CMakeLists.txt b/dumux/porousmediumflow/richards/CMakeLists.txt
index 1118c303d133c6316332ed8e359b45af1e556183..da07bb511a997fd78a7a1f16d8ab8fea7935f184 100644
--- a/dumux/porousmediumflow/richards/CMakeLists.txt
+++ b/dumux/porousmediumflow/richards/CMakeLists.txt
@@ -5,6 +5,5 @@ localresidual.hh
 model.hh
 newtonsolver.hh
 primaryvariableswitch.hh
-privarswitchnewtonsolver.hh
 volumevariables.hh
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/richards)
diff --git a/dumux/porousmediumflow/richards/model.hh b/dumux/porousmediumflow/richards/model.hh
index 57809e94d337726c62951cc4c5dfa74a4356ce58..a4eb5b344798472c2515b1ddd86eded09924d2c6 100644
--- a/dumux/porousmediumflow/richards/model.hh
+++ b/dumux/porousmediumflow/richards/model.hh
@@ -113,7 +113,6 @@
 #include "volumevariables.hh"
 #include "iofields.hh"
 #include "localresidual.hh"
-#include "primaryvariableswitch.hh"
 
 namespace Dumux {
 
@@ -239,10 +238,6 @@ public:
     using type = SwitchablePrimaryVariables<PrimaryVariablesVector, int>;
 };
 
-//! The primary variable switch for the richards model
-template<class TypeTag>
-struct PrimaryVariableSwitch<TypeTag, TTag::Richards> { using type = ExtendedRichardsPrimaryVariableSwitch; };
-
 /*!
  *\brief The fluid system used by the model.
  *
diff --git a/dumux/porousmediumflow/richards/newtonsolver.hh b/dumux/porousmediumflow/richards/newtonsolver.hh
index 60e6124f720a85e6c4a20bbd1295fe199269296c..cf69d3e025544698d9b482cfe59b13322ae3ec6f 100644
--- a/dumux/porousmediumflow/richards/newtonsolver.hh
+++ b/dumux/porousmediumflow/richards/newtonsolver.hh
@@ -25,7 +25,7 @@
 #define DUMUX_RICHARDS_NEWTON_SOLVER_HH
 
 #include <dumux/common/properties.hh>
-#include <dumux/porousmediumflow/richards/privarswitchnewtonsolver.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
 #include <dumux/discretization/elementsolution.hh>
 
 namespace Dumux {
@@ -39,15 +39,15 @@ namespace Dumux {
  * \todo make this typetag independent by extracting anything model specific from assembler
  *       or from possible ModelTraits.
  */
-template <class TypeTag, class Assembler, class LinearSolver>
-class RichardsNewtonSolver : public RichardsPrivarSwitchNewtonSolver<TypeTag, Assembler, LinearSolver>
+template <class Assembler, class LinearSolver>
+class RichardsNewtonSolver : public NewtonSolver<Assembler, LinearSolver>
 {
     using Scalar = typename Assembler::Scalar;
-    using ParentType = RichardsPrivarSwitchNewtonSolver<TypeTag, Assembler, LinearSolver>;
+    using ParentType = NewtonSolver<Assembler, LinearSolver>;
     using SolutionVector = typename Assembler::ResidualType;
 
     using MaterialLaw = typename Assembler::Problem::SpatialParams::MaterialLaw;
-    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+    using Indices = typename Assembler::GridVariables::VolumeVariables::Indices;
     enum { pressureIdx = Indices::pressureIdx };
 
 public:
diff --git a/dumux/porousmediumflow/richards/privarswitchnewtonsolver.hh b/dumux/porousmediumflow/richards/privarswitchnewtonsolver.hh
deleted file mode 100644
index 975273273c1a68460c53b2d7b3140d015f310c26..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/richards/privarswitchnewtonsolver.hh
+++ /dev/null
@@ -1,69 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- *   See the file COPYING for full copying permissions.                      *
- *                                                                           *
- *   This program is free software: you can redistribute it and/or modify    *
- *   it under the terms of the GNU General Public License as published by    *
- *   the Free Software Foundation, either version 2 of the License, or       *
- *   (at your option) any later version.                                     *
- *                                                                           *
- *   This program is distributed in the hope that it will be useful,         *
- *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
- *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
- *   GNU General Public License for more details.                            *
- *                                                                           *
- *   You should have received a copy of the GNU General Public License       *
- *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
- *****************************************************************************/
-/*!
- * \file
- * \ingroup RichardsModel
- * \brief A Richards model newton solver.
- */
-#ifndef DUMUX_RICHARDS_PRIVAR_SWITCH_NEWTON_SOLVER_HH
-#define DUMUX_RICHARDS_PRIVAR_SWITCH_NEWTON_SOLVER_HH
-
-#include <dumux/common/properties.hh>
-#include <dumux/nonlinear/newtonsolver.hh>
-#include <dumux/nonlinear/privarswitchnewtonsolver.hh>
-#include <dumux/discretization/elementsolution.hh>
-
-namespace Dumux {
-
-template <class TypeTag, class Assembler, class LinearSolver, bool enableWaterDiffusionInAir>
-class RichardsPrivarSwitchNewtonSolverImplementation;
-/*!
- * \ingroup RichardsModel
- * \brief A base for the richards newton solver which derives from the right base newton solver.
-  */
-template <class TypeTag, class Assembler, class LinearSolver>
-using RichardsPrivarSwitchNewtonSolver = RichardsPrivarSwitchNewtonSolverImplementation <TypeTag, Assembler, LinearSolver, getPropValue<TypeTag, Properties::EnableWaterDiffusionInAir>()>;
-
-/*!
- * \ingroup RichardsModel
- * \brief the case without a primary variables switch
-  */
-template <class TypeTag, class Assembler, class LinearSolver>
-class RichardsPrivarSwitchNewtonSolverImplementation<TypeTag, Assembler, LinearSolver, false> : public NewtonSolver<Assembler, LinearSolver>
-{
-    using ParentType = NewtonSolver<Assembler, LinearSolver>;
-public:
-    using ParentType::ParentType;
-};
-/*!
- * \ingroup RichardsModel
- * \brief the case with switchable primary variables
-  */
-template <class TypeTag, class Assembler, class LinearSolver>
-class RichardsPrivarSwitchNewtonSolverImplementation<TypeTag, Assembler, LinearSolver, true> : public PriVarSwitchNewtonSolver<Assembler, LinearSolver, GetPropType<TypeTag, Properties::PrimaryVariableSwitch>>
-{
-    using PrimaryVariableSwitch = GetPropType<TypeTag, Properties::PrimaryVariableSwitch>;
-    using ParentType = PriVarSwitchNewtonSolver<Assembler, LinearSolver, PrimaryVariableSwitch>;
-public:
-    using ParentType::ParentType;
-};
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/porousmediumflow/richards/volumevariables.hh b/dumux/porousmediumflow/richards/volumevariables.hh
index fe83f039ccd219eaf922929b57187d512759c689..3b629d328af9c29ec7eddb82b945d6d9d8efea13 100644
--- a/dumux/porousmediumflow/richards/volumevariables.hh
+++ b/dumux/porousmediumflow/richards/volumevariables.hh
@@ -34,8 +34,21 @@
 #include <dumux/material/constants.hh>
 #include <dumux/material/solidstates/updatesolidvolumefractions.hh>
 
+#include "primaryvariableswitch.hh"
+
 namespace Dumux {
 
+namespace Detail {
+//! helper structs to conditionally use a primary variable switch or not
+struct VolVarsWithPVSwitch
+{
+    using PrimaryVariableSwitch = ExtendedRichardsPrimaryVariableSwitch;
+};
+
+struct VolVarsWithOutPVSwitch
+{};
+}
+
 /*!
  * \ingroup RichardsModel
  * \brief Volume averaged quantities required by the Richards model.
@@ -47,6 +60,8 @@ template <class Traits>
 class RichardsVolumeVariables
 : public PorousMediumFlowVolumeVariables<Traits>
 , public EnergyVolumeVariables<Traits, RichardsVolumeVariables<Traits> >
+, public std::conditional_t<Traits::ModelTraits::enableMolecularDiffusion(),
+                            Detail::VolVarsWithPVSwitch, Detail::VolVarsWithOutPVSwitch>
 {
     using ParentType = PorousMediumFlowVolumeVariables<Traits>;
     using EnergyVolVars = EnergyVolumeVariables<Traits, RichardsVolumeVariables<Traits> >;
diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/main.cc b/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/main.cc
index c4b80552e9dd69f49956330bfc77a8e62a96ad92..5024433d0c2522f5a1603618d86110037bb5a335 100644
--- a/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/main.cc
+++ b/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/main.cc
@@ -43,7 +43,7 @@
 
 #include <dumux/multidomain/staggeredtraits.hh>
 #include <dumux/multidomain/fvassembler.hh>
-#include <dumux/multidomain/privarswitchnewtonsolver.hh>
+#include <dumux/multidomain/newtonsolver.hh>
 
 #include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
 
@@ -202,11 +202,8 @@ int main(int argc, char** argv) try
     using LinearSolver = UMFPackBackend;
     auto linearSolver = std::make_shared<LinearSolver>();
 
-    // the primary variable switches used by the sub models
-    using PriVarSwitchTuple = std::tuple<NoPrimaryVariableSwitch, NoPrimaryVariableSwitch, GetPropType<DarcyTypeTag, Properties::PrimaryVariableSwitch>>;
-
     // the non-linear solver
-    using NewtonSolver = MultiDomainPriVarSwitchNewtonSolver<Assembler, LinearSolver, CouplingManager, PriVarSwitchTuple>;
+    using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
     NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
 
     // time loop
diff --git a/test/porousmediumflow/2p1c/implicit/main.cc b/test/porousmediumflow/2p1c/implicit/main.cc
index 5d0069b9ea9c469d357513255c015d3fa1b32db3..53d3209a6c444a596b13ef10f9bf2a40072ffa20 100644
--- a/test/porousmediumflow/2p1c/implicit/main.cc
+++ b/test/porousmediumflow/2p1c/implicit/main.cc
@@ -42,7 +42,7 @@
 #include <dumux/common/defaultusagemessage.hh>
 
 #include <dumux/linear/seqsolverbackend.hh>
-#include <dumux/nonlinear/privarswitchnewtonsolver.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
 
 #include <dumux/assembly/fvassembler.hh>
 
@@ -128,8 +128,7 @@ int main(int argc, char** argv) try
     auto linearSolver = std::make_shared<LinearSolver>();
 
     // the non-linear solver
-    using PrimaryVariableSwitch = GetPropType<TypeTag, Properties::PrimaryVariableSwitch>;
-    using NewtonSolver = Dumux::PriVarSwitchNewtonSolver<Assembler, LinearSolver, PrimaryVariableSwitch>;
+    using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
 
     // time loop
diff --git a/test/porousmediumflow/2p2c/implicit/injection/main.cc b/test/porousmediumflow/2p2c/implicit/injection/main.cc
index b280d8e7e43315b9cafe74257a5da5eb81acd546..22041b741c4c5287a63bec961e94191250b267fa 100644
--- a/test/porousmediumflow/2p2c/implicit/injection/main.cc
+++ b/test/porousmediumflow/2p2c/implicit/injection/main.cc
@@ -38,7 +38,7 @@
 #include <dumux/common/defaultusagemessage.hh>
 
 #include <dumux/linear/amgbackend.hh>
-#include <dumux/nonlinear/privarswitchnewtonsolver.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
 
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/assembly/diffmethod.hh>
@@ -141,8 +141,7 @@ int main(int argc, char** argv) try
     auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
 
     // the non-linear solver
-    using NewtonSolver = PriVarSwitchNewtonSolver<Assembler, LinearSolver,
-                                                 GetPropType<TypeTag, Properties::PrimaryVariableSwitch>>;
+    using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
 
     // time loop
diff --git a/test/porousmediumflow/2p2c/implicit/mpnccomparison/main.cc b/test/porousmediumflow/2p2c/implicit/mpnccomparison/main.cc
index d67cd035fd7477654d1e8792ef45ed412cc14cdf..cbc85db68595b8c2149c882acdd32db917da471b 100644
--- a/test/porousmediumflow/2p2c/implicit/mpnccomparison/main.cc
+++ b/test/porousmediumflow/2p2c/implicit/mpnccomparison/main.cc
@@ -38,7 +38,7 @@
 #include <dumux/common/defaultusagemessage.hh>
 
 #include <dumux/linear/amgbackend.hh>
-#include <dumux/nonlinear/privarswitchnewtonsolver.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
 
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/assembly/diffmethod.hh>
@@ -152,8 +152,7 @@ int main(int argc, char** argv) try
     auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
 
     // the non-linear solver
-    using NewtonSolver = PriVarSwitchNewtonSolver<Assembler, LinearSolver,
-                                                 GetPropType<TypeTag, Properties::PrimaryVariableSwitch>>;
+    using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
 
     // time loop
diff --git a/test/porousmediumflow/2p2c/implicit/waterair/main.cc b/test/porousmediumflow/2p2c/implicit/waterair/main.cc
index 6a6e60ebae5c5c854d30b39847ff4930af95c18b..bd87600a02a7ed530e385c62d8439ef2c5781ce3 100644
--- a/test/porousmediumflow/2p2c/implicit/waterair/main.cc
+++ b/test/porousmediumflow/2p2c/implicit/waterair/main.cc
@@ -38,7 +38,7 @@
 #include <dumux/common/defaultusagemessage.hh>
 
 #include <dumux/linear/amgbackend.hh>
-#include <dumux/nonlinear/privarswitchnewtonsolver.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
 
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/assembly/diffmethod.hh>
@@ -141,8 +141,7 @@ int main(int argc, char** argv) try
     auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
 
     // the non-linear solver
-    using NewtonSolver = PriVarSwitchNewtonSolver<Assembler, LinearSolver,
-                                                 GetPropType<TypeTag, Properties::PrimaryVariableSwitch>>;
+    using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
 
     // time loop
diff --git a/test/porousmediumflow/2pnc/implicit/diffusion/main.cc b/test/porousmediumflow/2pnc/implicit/diffusion/main.cc
index ce214776ebed8b7d96dc27b0d1c4b3ced2927a5f..d5bd4a9e7fef321ce95eb62f9309fda1fa73a974 100644
--- a/test/porousmediumflow/2pnc/implicit/diffusion/main.cc
+++ b/test/porousmediumflow/2pnc/implicit/diffusion/main.cc
@@ -41,7 +41,7 @@
 #include <dumux/common/defaultusagemessage.hh>
 
 #include <dumux/linear/amgbackend.hh>
-#include <dumux/nonlinear/privarswitchnewtonsolver.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
 
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/assembly/diffmethod.hh>
@@ -149,8 +149,7 @@ int main(int argc, char** argv) try
     auto linearSolver =  std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
 
     // the non-linear solver
-    using NewtonSolver = PriVarSwitchNewtonSolver<Assembler, LinearSolver,
-                                                  GetPropType<TypeTag, Properties::PrimaryVariableSwitch>>;
+    using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
 
     // time loop
diff --git a/test/porousmediumflow/2pnc/implicit/fuelcell/main.cc b/test/porousmediumflow/2pnc/implicit/fuelcell/main.cc
index d2baee027e73c14e3757abac64b36b0bffc53ac3..be3c5f600f656ed2c3f0fc77eb1b900866264467 100644
--- a/test/porousmediumflow/2pnc/implicit/fuelcell/main.cc
+++ b/test/porousmediumflow/2pnc/implicit/fuelcell/main.cc
@@ -39,7 +39,7 @@
 #include <dumux/common/defaultusagemessage.hh>
 
 #include <dumux/linear/amgbackend.hh>
-#include <dumux/nonlinear/privarswitchnewtonsolver.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
 
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/assembly/diffmethod.hh>
@@ -150,8 +150,7 @@ int main(int argc, char** argv) try
     auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
 
     // the non-linear solver
-    using NewtonSolver = PriVarSwitchNewtonSolver<Assembler, LinearSolver,
-                                                  GetPropType<TypeTag, Properties::PrimaryVariableSwitch>>;
+    using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
 
     // time loop
diff --git a/test/porousmediumflow/2pncmin/implicit/main.cc b/test/porousmediumflow/2pncmin/implicit/main.cc
index 49b82846aa74d42b1950c7e87d57cbcf63fb0fed..6f2f1e972448d2fff7182d1cdf5a7df34dafb02d 100644
--- a/test/porousmediumflow/2pncmin/implicit/main.cc
+++ b/test/porousmediumflow/2pncmin/implicit/main.cc
@@ -41,7 +41,7 @@
 #include <dumux/common/defaultusagemessage.hh>
 
 #include <dumux/linear/amgbackend.hh>
-#include <dumux/nonlinear/privarswitchnewtonsolver.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
 
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/assembly/diffmethod.hh>
@@ -169,8 +169,7 @@ int main(int argc, char** argv) try
     auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
 
     // the non-linear solver
-    using NewtonSolver = PriVarSwitchNewtonSolver<Assembler, LinearSolver,
-                                                  GetPropType<TypeTag, Properties::PrimaryVariableSwitch>>;
+    using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
 
     // time loop
diff --git a/test/porousmediumflow/3p3c/implicit/columnxylol/main.cc b/test/porousmediumflow/3p3c/implicit/columnxylol/main.cc
index 9bf46173b7c0fbd238e7f877e19cc7b3d6c19f02..062be07df41594e70fc10ba2d9736817a781acb7 100644
--- a/test/porousmediumflow/3p3c/implicit/columnxylol/main.cc
+++ b/test/porousmediumflow/3p3c/implicit/columnxylol/main.cc
@@ -37,7 +37,7 @@
 #include <dumux/common/timeloop.hh>
 
 #include <dumux/linear/amgbackend.hh>
-#include <dumux/nonlinear/privarswitchnewtonsolver.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
 
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/assembly/diffmethod.hh>
@@ -155,8 +155,7 @@ int main(int argc, char** argv) try
     auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
 
     // the non-linear solver
-    using NewtonSolver = PriVarSwitchNewtonSolver<Assembler, LinearSolver,
-                                                  GetPropType<TypeTag, Properties::PrimaryVariableSwitch>>;
+    using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
 
     // time loop
diff --git a/test/porousmediumflow/3p3c/implicit/infiltration/main.cc b/test/porousmediumflow/3p3c/implicit/infiltration/main.cc
index 9bf46173b7c0fbd238e7f877e19cc7b3d6c19f02..062be07df41594e70fc10ba2d9736817a781acb7 100644
--- a/test/porousmediumflow/3p3c/implicit/infiltration/main.cc
+++ b/test/porousmediumflow/3p3c/implicit/infiltration/main.cc
@@ -37,7 +37,7 @@
 #include <dumux/common/timeloop.hh>
 
 #include <dumux/linear/amgbackend.hh>
-#include <dumux/nonlinear/privarswitchnewtonsolver.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
 
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/assembly/diffmethod.hh>
@@ -155,8 +155,7 @@ int main(int argc, char** argv) try
     auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
 
     // the non-linear solver
-    using NewtonSolver = PriVarSwitchNewtonSolver<Assembler, LinearSolver,
-                                                  GetPropType<TypeTag, Properties::PrimaryVariableSwitch>>;
+    using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
 
     // time loop
diff --git a/test/porousmediumflow/3p3c/implicit/kuevette/main.cc b/test/porousmediumflow/3p3c/implicit/kuevette/main.cc
index 9bf46173b7c0fbd238e7f877e19cc7b3d6c19f02..062be07df41594e70fc10ba2d9736817a781acb7 100644
--- a/test/porousmediumflow/3p3c/implicit/kuevette/main.cc
+++ b/test/porousmediumflow/3p3c/implicit/kuevette/main.cc
@@ -37,7 +37,7 @@
 #include <dumux/common/timeloop.hh>
 
 #include <dumux/linear/amgbackend.hh>
-#include <dumux/nonlinear/privarswitchnewtonsolver.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
 
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/assembly/diffmethod.hh>
@@ -155,8 +155,7 @@ int main(int argc, char** argv) try
     auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
 
     // the non-linear solver
-    using NewtonSolver = PriVarSwitchNewtonSolver<Assembler, LinearSolver,
-                                                  GetPropType<TypeTag, Properties::PrimaryVariableSwitch>>;
+    using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
 
     // time loop
diff --git a/test/porousmediumflow/3pwateroil/implicit/main.cc b/test/porousmediumflow/3pwateroil/implicit/main.cc
index 5a50679a5462cbf39eb06c8f62ca20b1fdbef1d1..475f74fde915078a72de18a837aef772a3dfb0b8 100644
--- a/test/porousmediumflow/3pwateroil/implicit/main.cc
+++ b/test/porousmediumflow/3pwateroil/implicit/main.cc
@@ -37,7 +37,7 @@
 #include <dumux/common/timeloop.hh>
 
 #include <dumux/linear/amgbackend.hh>
-#include <dumux/nonlinear/privarswitchnewtonsolver.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
 
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/assembly/diffmethod.hh>
@@ -154,8 +154,7 @@ int main(int argc, char** argv) try
     auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
 
     // the non-linear solver
-    using NewtonSolver = PriVarSwitchNewtonSolver<Assembler, LinearSolver,
-                                                  GetPropType<TypeTag, Properties::PrimaryVariableSwitch>>;
+    using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
 
     // time loop
diff --git a/test/porousmediumflow/co2/implicit/main.cc b/test/porousmediumflow/co2/implicit/main.cc
index 7a51ffd9e5a8f491dbd5d5a9ece3d713434d5ec6..fa9ee240977fb045664051e07d6de59347f76a0f 100644
--- a/test/porousmediumflow/co2/implicit/main.cc
+++ b/test/porousmediumflow/co2/implicit/main.cc
@@ -38,7 +38,7 @@
 #include <dumux/common/defaultusagemessage.hh>
 
 #include <dumux/linear/amgbackend.hh>
-#include <dumux/nonlinear/privarswitchnewtonsolver.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
 
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/assembly/diffmethod.hh>
@@ -131,8 +131,7 @@ int main(int argc, char** argv) try
     auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
 
     // the non-linear solver
-    using NewtonSolver = PriVarSwitchNewtonSolver<Assembler, LinearSolver,
-                                                  GetPropType<TypeTag, Properties::PrimaryVariableSwitch>>;
+    using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
 
     // time loop
diff --git a/test/porousmediumflow/richards/implicit/analytical/main.cc b/test/porousmediumflow/richards/implicit/analytical/main.cc
index dc266ba1917fe703c4002c4c49ff27b34e21b2e8..bd09674a53879f9a56d4212e07677cdbc6ca25cf 100644
--- a/test/porousmediumflow/richards/implicit/analytical/main.cc
+++ b/test/porousmediumflow/richards/implicit/analytical/main.cc
@@ -151,7 +151,7 @@ int main(int argc, char** argv) try
     auto linearSolver = std::make_shared<LinearSolver>();
 
     // the non-linear solver
-    using NewtonSolver = Dumux::RichardsNewtonSolver<TypeTag, Assembler, LinearSolver>;
+    using NewtonSolver = Dumux::RichardsNewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
 
     // time loop
diff --git a/test/porousmediumflow/richards/implicit/lens/main.cc b/test/porousmediumflow/richards/implicit/lens/main.cc
index 06454b52d8de1c027d0332e0679b00a4b11e55eb..dd88b9298def3e0eafe3c41dbab8be6e708036d0 100644
--- a/test/porousmediumflow/richards/implicit/lens/main.cc
+++ b/test/porousmediumflow/richards/implicit/lens/main.cc
@@ -143,7 +143,7 @@ int main(int argc, char** argv) try
     auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
 
     // the non-linear solver
-    using NewtonSolver = Dumux::RichardsNewtonSolver<TypeTag, Assembler, LinearSolver>;
+    using NewtonSolver = Dumux::RichardsNewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
 
     // time loop
diff --git a/test/porousmediumflow/richards/implicit/nonisothermal/conduction/main.cc b/test/porousmediumflow/richards/implicit/nonisothermal/conduction/main.cc
index 0123cf6476209960632f02a27321934e3a0aaabb..9afd079379eb7f8a9c88a890faaf10335b29dcf3 100644
--- a/test/porousmediumflow/richards/implicit/nonisothermal/conduction/main.cc
+++ b/test/porousmediumflow/richards/implicit/nonisothermal/conduction/main.cc
@@ -152,7 +152,7 @@ int main(int argc, char** argv) try
     auto linearSolver = std::make_shared<LinearSolver>();
 
     // the non-linear solver
-    using NewtonSolver = Dumux::RichardsNewtonSolver<TypeTag, Assembler, LinearSolver>;
+    using NewtonSolver = Dumux::RichardsNewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
 
     // time loop
diff --git a/test/porousmediumflow/richards/implicit/nonisothermal/convection/main.cc b/test/porousmediumflow/richards/implicit/nonisothermal/convection/main.cc
index 0123cf6476209960632f02a27321934e3a0aaabb..9afd079379eb7f8a9c88a890faaf10335b29dcf3 100644
--- a/test/porousmediumflow/richards/implicit/nonisothermal/convection/main.cc
+++ b/test/porousmediumflow/richards/implicit/nonisothermal/convection/main.cc
@@ -152,7 +152,7 @@ int main(int argc, char** argv) try
     auto linearSolver = std::make_shared<LinearSolver>();
 
     // the non-linear solver
-    using NewtonSolver = Dumux::RichardsNewtonSolver<TypeTag, Assembler, LinearSolver>;
+    using NewtonSolver = Dumux::RichardsNewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
 
     // time loop
diff --git a/test/porousmediumflow/richards/implicit/nonisothermal/evaporation/main.cc b/test/porousmediumflow/richards/implicit/nonisothermal/evaporation/main.cc
index 37f264d262047ec5932cfa21f553499814c79aa4..720d855c0199debaf53a650d0cede0c56b063557 100644
--- a/test/porousmediumflow/richards/implicit/nonisothermal/evaporation/main.cc
+++ b/test/porousmediumflow/richards/implicit/nonisothermal/evaporation/main.cc
@@ -124,7 +124,7 @@ int main(int argc, char** argv) try
     auto linearSolver = std::make_shared<LinearSolver>();
 
     // the non-linear solver
-    using NewtonSolver = Dumux::RichardsNewtonSolver<TypeTag, Assembler, LinearSolver>;
+    using NewtonSolver = RichardsNewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
 
     // time loop
diff --git a/test/porousmediumflow/richardsnc/implicit/main.cc b/test/porousmediumflow/richardsnc/implicit/main.cc
index 02ba9a34b889261a05d90d2dcb09c2fa5f2b80db..8f278f70db9db178d3e17dd94be52d7d86cad2a3 100644
--- a/test/porousmediumflow/richardsnc/implicit/main.cc
+++ b/test/porousmediumflow/richardsnc/implicit/main.cc
@@ -152,7 +152,7 @@ int main(int argc, char** argv) try
     auto linearSolver = std::make_shared<LinearSolver>();
 
     // the non-linear solver
-    using NewtonSolver = Dumux::RichardsNewtonSolver<TypeTag, Assembler, LinearSolver>;
+    using NewtonSolver = Dumux::RichardsNewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
 
     // time loop