diff --git a/dumux/discretization/extrusion.hh b/dumux/discretization/extrusion.hh
index 758c76d8c26a96782e53c874d2d44e7e90b864ff..8832fa5fef4987e5277da8e5d426dee5afb441e8 100644
--- a/dumux/discretization/extrusion.hh
+++ b/dumux/discretization/extrusion.hh
@@ -53,10 +53,13 @@ struct NoExtrusion
 /*!
  * \ingroup Discretization
  * \brief Rotation symmetric extrusion policy for rotating about an external axis
+ * \tparam radAx The radial axis perpendicular to the symmetry axis (0 = x, 1 = y)
  */
-template<int radialAxis = 0>
+template<int radAx = 0>
 struct RotationalExtrusion
 {
+    static constexpr int radialAxis = radAx;
+
     /*!
      * \brief Transformed sub-control-volume face area
      * \note Mid-point rule integrals are only exact for constants
diff --git a/dumux/freeflow/navierstokes/staggered/localresidual.hh b/dumux/freeflow/navierstokes/staggered/localresidual.hh
index 64083aa4584dca6e024087ff2c6e534896262ecf..6baa841272ba21006e8d85be608344bf15cfe85d 100644
--- a/dumux/freeflow/navierstokes/staggered/localresidual.hh
+++ b/dumux/freeflow/navierstokes/staggered/localresidual.hh
@@ -169,6 +169,24 @@ public:
         source += problem.gravity()[scvf.directionIndex()] * insideVolVars.density();
         source += problem.source(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[Indices::velocity(scvf.directionIndex())];
 
+        // Axisymmetric problems in 2D feature an extra source terms arising from the transformation to cylindrical coordinates.
+        // See Ferziger/Peric: Computational methods for fluid dynamics chapter 8.
+        // https://doi.org/10.1007/978-3-540-68228-8 (page 301)
+        if constexpr (ModelTraits::dim() == 2 && std::is_same_v<Extrusion, RotationalExtrusion<Extrusion::radialAxis>>)
+        {
+            if (scvf.directionIndex() == Extrusion::radialAxis)
+            {
+                // Velocity term
+                const auto r = scvf.center()[scvf.directionIndex()];
+                source -= -2.0*insideVolVars.effectiveViscosity() * elemFaceVars[scvf].velocitySelf() / (r*r);
+
+                // Pressure term (needed because we incorporate pressure in terms of a surface integral).
+                // grad(p) becomes div(pI) + (p/r)*n_r in cylindrical coordinates. The second term
+                // is new with respect to Cartesian coordinates and handled below as a source term.
+                source += insideVolVars.pressure()/r;
+            }
+        }
+
         return source;
     }
 
diff --git a/test/freeflow/navierstokes/channel/pipe/main.cc b/test/freeflow/navierstokes/channel/pipe/main.cc
index 4ca8e828d689d6f4dbf6c4ed70ca03b0f671fa86..4f9691791a31279310161883e4e4f64e23219ae7 100644
--- a/test/freeflow/navierstokes/channel/pipe/main.cc
+++ b/test/freeflow/navierstokes/channel/pipe/main.cc
@@ -101,7 +101,6 @@ int main(int argc, char** argv) try
     SolutionVector sol;
     sol[GridGeometry::cellCenterIdx()].resize(gridGeometry->numCellCenterDofs());
     sol[GridGeometry::faceIdx()].resize(gridGeometry->numFaceDofs());
-    problem->applyInitialSolution(sol);
 
     // the grid variables
     using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
diff --git a/test/freeflow/navierstokes/channel/pipe/params.input b/test/freeflow/navierstokes/channel/pipe/params.input
index 268e5b9bf993d9f7749869b1643baabe47f31725..9bcb4dacc3c10c43bbeba8a960d68da7c40292c3 100644
--- a/test/freeflow/navierstokes/channel/pipe/params.input
+++ b/test/freeflow/navierstokes/channel/pipe/params.input
@@ -20,3 +20,6 @@ LiquidDensity = 1000
 
 [Vtk]
 AddVelocity = 1
+
+[Assembly]
+NumericDifference.BaseEpsilon =  1e-3
diff --git a/test/freeflow/navierstokes/channel/pipe/problem.hh b/test/freeflow/navierstokes/channel/pipe/problem.hh
index cf5cc443b91072fbd58f8f1248d05b36e468708f..91668701cdd87676fa18c63501cc941ccc36c94d 100644
--- a/test/freeflow/navierstokes/channel/pipe/problem.hh
+++ b/test/freeflow/navierstokes/channel/pipe/problem.hh
@@ -109,7 +109,7 @@ public:
     { return analyticalSolution(globalPos); }
 
     PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
-    { return analyticalSolution(globalPos); }
+    { return PrimaryVariables(0.0); }
 
     PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const
     {
@@ -120,7 +120,7 @@ public:
         const auto y = globalPos[1] - this->gridGeometry().bBoxMin()[1];
         values[Indices::velocityXIdx] = 0.0;
         values[Indices::velocityYIdx] = 2.0*meanInletVelocity_*(1.0 - r*r/(pipeRadius_*pipeRadius_));
-        values[Indices::pressureIdx] = 1e5 + (pipeLength_-y)*meanInletVelocity_*8.0*mu_/(pipeRadius_*pipeRadius_);
+        values[Indices::pressureIdx] = (pipeLength_-y)*meanInletVelocity_*8.0*mu_/(pipeRadius_*pipeRadius_);
 
         return values;
     }