diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh
index b33f113a46d5ad5f16f7c117120f02b620dfeda5..647a0c6ba28ff9b7ad92e0d969a1002bab4d9b8f 100644
--- a/dumux/nonlinear/newtonsolver.hh
+++ b/dumux/nonlinear/newtonsolver.hh
@@ -444,12 +444,18 @@ public:
         {
             if (numSteps_ == 0)
             {
-                Scalar norm2 = b.two_norm2();
-                if (comm_.size() > 1)
-                    norm2 = comm_.sum(norm2);
+                if constexpr (Detail::hasNorm<LinearSolver, SolutionVector>())
+                    initialResidual_ = this->linearSolver().norm(b);
 
-                using std::sqrt;
-                initialResidual_ = sqrt(norm2);
+                else
+                {
+                    Scalar norm2 = b.two_norm2();
+                    if (comm_.size() > 1)
+                        norm2 = comm_.sum(norm2);
+
+                    using std::sqrt;
+                    initialResidual_ = sqrt(norm2);
+                }
             }
 
             // solve by calling the appropriate implementation depending on whether the linear solver
@@ -789,7 +795,10 @@ protected:
     void computeResidualReduction_(const SolutionVector &uCurrentIter)
     {
         if constexpr (Detail::hasNorm<LinearSolver, SolutionVector>())
-            residualNorm_ = this->linearSolver().norm(uCurrentIter);
+        {
+            this->assembler().assembleResidual(uCurrentIter);
+            residualNorm_ = this->linearSolver().norm(this->assembler().residual());
+        }
         else
             residualNorm_ = this->assembler().residualNorm(uCurrentIter);