diff --git a/test/freeflow/navierstokes/permeabilityupscaling/main.cc b/test/freeflow/navierstokes/permeabilityupscaling/main.cc
index b80a3c7e9706b671cdbe0787b1a7d69651aac7c1..3480a1396a602de32783da5cdedb76d9fbadf8fa 100644
--- a/test/freeflow/navierstokes/permeabilityupscaling/main.cc
+++ b/test/freeflow/navierstokes/permeabilityupscaling/main.cc
@@ -19,7 +19,6 @@
 #include <dumux/common/initialize.hh>
 #include <dumux/common/dumuxmessage.hh>
 #include <dumux/common/parameters.hh>
-#include <dumux/common/indextraits.hh>
 #include <dumux/common/properties.hh>
 
 #include <dumux/io/vtkoutputmodule.hh>
@@ -27,8 +26,7 @@
 #include <dumux/io/grid/gridmanager_yasp.hh>
 
 #include <dumux/linear/istlsolvers.hh>
-#include <dumux/linear/linearsolvertraits.hh>
-#include <dumux/linear/linearalgebratraits.hh>
+#include <dumux/linear/stokes_solver.hh>
 
 #include <dumux/multidomain/fvassembler.hh>
 #include <dumux/multidomain/traits.hh>
@@ -64,6 +62,39 @@ auto makeAveragingVolume(const Dune::FieldVector<ctype, dim>& center,
     };
 }
 
+template<class Vector, class MomGG, class MassGG, class MomP, class MomIdx, class MassIdx>
+auto dirichletDofs(std::shared_ptr<MomGG> momentumGridGeometry,
+                   std::shared_ptr<MassGG> massGridGeometry,
+                   std::shared_ptr<MomP> momentumProblem,
+                   MomIdx momentumIdx, MassIdx massIdx)
+{
+    Vector dirichletDofs;
+    dirichletDofs[momentumIdx].resize(momentumGridGeometry->numDofs());
+    dirichletDofs[massIdx].resize(massGridGeometry->numDofs());
+    dirichletDofs = 0.0;
+
+    auto fvGeometry = localView(*momentumGridGeometry);
+    for (const auto& element : elements(momentumGridGeometry->gridView()))
+    {
+        fvGeometry.bind(element);
+        for (const auto& scvf : scvfs(fvGeometry))
+        {
+            if (!scvf.boundary() || !scvf.isFrontal())
+                continue;
+            const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
+            if (scv.boundary())
+            {
+                const auto bcTypes = momentumProblem->boundaryTypes(element, scvf);
+                for (int i = 0; i < bcTypes.size(); ++i)
+                    if (bcTypes.isDirichlet(i))
+                        dirichletDofs[momentumIdx][scv.dofIndex()][i] = 1.0;
+            }
+        }
+    }
+
+    return dirichletDofs;
+}
+
 int main(int argc, char** argv)
 {
     using namespace Dumux;
@@ -142,14 +173,26 @@ int main(int argc, char** argv)
     vtkWriter.addVelocityOutput(std::make_shared<NavierStokesVelocityOutput<MassGridVariables>>());
     vtkWriter.write(0.0);
 
-    // the linear solver
-    using LinearSolver = UMFPackIstlSolver<SeqLinearSolverTraits, LinearAlgebraTraitsFromAssembler<Assembler>>;
-    auto linearSolver = std::make_shared<LinearSolver>();
-
-    // solve the non-linear system
-    using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
-    NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
-    nonLinearSolver.solve(x);
+    // the linearize and solve
+    if (getParam<bool>("LinearSolver.UseIterativeSolver", true))
+    {
+        using Matrix = typename Assembler::JacobianMatrix;
+        using Vector = typename Assembler::ResidualType;
+        using LinearSolver = StokesSolver<Matrix, Vector, MomentumGridGeometry, MassGridGeometry>;
+        auto dDofs = dirichletDofs<Vector>(momentumGridGeometry, massGridGeometry, momentumProblem, momentumIdx, massIdx);
+        auto linearSolver = std::make_shared<LinearSolver>(momentumGridGeometry, massGridGeometry, dDofs);
+        using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
+        NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+        nonLinearSolver.solve(x);
+    }
+    else
+    {
+        using LinearSolver = UMFPackIstlSolver<SeqLinearSolverTraits, LinearAlgebraTraitsFromAssembler<Assembler>>;
+        auto linearSolver = std::make_shared<LinearSolver>();
+        using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
+        NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+        nonLinearSolver.solve(x);
+    }
 
     // post-processing and output
     vtkWriter.write(1.0);
diff --git a/test/freeflow/navierstokes/permeabilityupscaling/params.input b/test/freeflow/navierstokes/permeabilityupscaling/params.input
index 24b2892cfccc063fbaab3f654b8d04b0b3f379e9..770c893fe575770e5a88df881d808bc8f9e58cde 100644
--- a/test/freeflow/navierstokes/permeabilityupscaling/params.input
+++ b/test/freeflow/navierstokes/permeabilityupscaling/params.input
@@ -11,12 +11,21 @@ EnableInertiaTerms = false
 [Component]
 # this test requires a unit fluid
 LiquidDensity = 1
-LiquidKinematicViscosity = 1
+LiquidDynamicViscosity = 1
 
 [ Newton ]
 MaxSteps = 20
 MaxRelativeShift = 1e-8
 
+[LinearSolver]
+MaxIterations = 500
+ResidualReduction = 1e-10
+SymmetrizeDirichlet = true
+DirectSolverForVelocity = false
+GMResRestart = 500
+Type = gmres
+Verbosity = 1
+
 [Vtk]
 WriteFaceData = false