diff --git a/MathLib/LinAlg/LinAlg.h b/MathLib/LinAlg/LinAlg.h
index 0429c581cf4..0a587933c66 100644
--- a/MathLib/LinAlg/LinAlg.h
+++ b/MathLib/LinAlg/LinAlg.h
@@ -11,6 +11,7 @@
 #pragma once
 
 #include <cassert>
+
 #include "BaseLib/Error.h"
 #include "LinAlgEnums.h"
 
@@ -269,3 +270,47 @@ void finalizeAssembly(EigenVector& A);
 } // namespace MathLib
 
 #endif
+
+namespace MathLib::LinAlg
+{
+
+/**
+ * Compute the relative norm between two vectors by \f$ \|x-y\|/\|x\| \f$.
+ *
+ * \param x Vector x
+ * \param y Vector y
+ * \param norm_type The norm type of global vector
+ * \return \f$ \|x-y\|/\|x\| \f$.
+ */
+template <typename VectorType>
+double computeRelativeNorm(VectorType const& x, VectorType const& y,
+                           MathLib::VecNormType norm_type)
+{
+    if (norm_type == MathLib::VecNormType::INVALID)
+    {
+        OGS_FATAL("An invalid norm type has been passed");
+    }
+
+    // Stores \f$diff = x-y\f$.
+    VectorType diff;
+    MathLib::LinAlg::copy(x, diff);
+    MathLib::LinAlg::axpy(diff, -1.0, y);
+
+    const double norm_diff = MathLib::LinAlg::norm(diff, norm_type);
+
+    const double norm_x = MathLib::LinAlg::norm(x, norm_type);
+    if (norm_x > std::numeric_limits<double>::epsilon())
+    {
+        return norm_diff / norm_x;
+    }
+
+    // Both of norm_x and norm_diff are close to zero
+    if (norm_diff < std::numeric_limits<double>::epsilon())
+    {
+        return 1.0;
+    }
+
+    // Only norm_x is close to zero
+    return norm_diff / std::numeric_limits<double>::epsilon();
+}
+}  // namespace MathLib::LinAlg
diff --git a/NumLib/ODESolver/TimeDiscretization.cpp b/NumLib/ODESolver/TimeDiscretization.cpp
index 9a7d414f873..5025ff4d6c2 100644
--- a/NumLib/ODESolver/TimeDiscretization.cpp
+++ b/NumLib/ODESolver/TimeDiscretization.cpp
@@ -15,40 +15,6 @@
 
 namespace NumLib
 {
-
-double computeRelativeNorm(GlobalVector const& x,
-                           GlobalVector const& y,
-                           MathLib::VecNormType norm_type)
-{
-    if (norm_type == MathLib::VecNormType::INVALID)
-    {
-        OGS_FATAL("An invalid norm type has been passed");
-    }
-
-    // Stores \f$ u_{n+1}-u_{n}\f$.
-    GlobalVector dx;
-    MathLib::LinAlg::copy(x, dx);  // copy x to dx.
-
-    // dx = x - y --> x - dx --> dx
-    MathLib::LinAlg::axpy(dx, -1.0, y);
-    const double norm_dx = MathLib::LinAlg::norm(dx, norm_type);
-
-    const double norm_x = MathLib::LinAlg::norm(x, norm_type);
-    if (norm_x > std::numeric_limits<double>::epsilon())
-    {
-        return norm_dx / norm_x;
-    }
-
-    // Both of norm_x and norm_dx are close to zero
-    if (norm_dx < std::numeric_limits<double>::epsilon())
-    {
-        return 1.0;
-    }
-
-    // Only norm_x is close to zero
-    return norm_dx / std::numeric_limits<double>::epsilon();
-}
-
 void BackwardEuler::getWeightedOldX(GlobalVector& y,
                                     GlobalVector const& x_old) const
 {
diff --git a/NumLib/ODESolver/TimeDiscretization.h b/NumLib/ODESolver/TimeDiscretization.h
index 4e1363e9c5c..d96ef4b058a 100644
--- a/NumLib/ODESolver/TimeDiscretization.h
+++ b/NumLib/ODESolver/TimeDiscretization.h
@@ -18,17 +18,6 @@
 
 namespace NumLib
 {
-/**
- * Compute the relative norm between two vectors by \f$ \|x-y\|/\|x\| \f$.
- *
- * \param x Vector x
- * \param y Vector y
- * \param norm_type The norm type of global vector
- * \return \f$ \|x-y\|/\|x\| \f$.
- */
-double computeRelativeNorm(GlobalVector const& x, GlobalVector const& y,
-                           MathLib::VecNormType norm_type);
-
 //! \addtogroup ODESolver
 //! @{
 
diff --git a/ProcessLib/TimeLoop.cpp b/ProcessLib/TimeLoop.cpp
index 3138f988f71..99675d1c775 100644
--- a/ProcessLib/TimeLoop.cpp
+++ b/ProcessLib/TimeLoop.cpp
@@ -332,7 +332,7 @@ std::pair<double, bool> TimeLoop::computeTimeStepping(
 
         const double solution_error =
             computationOfChangeNeeded(timestep_algorithm, t)
-                ? NumLib::computeRelativeNorm(
+                ? MathLib::LinAlg::computeRelativeNorm(
                       x, x_prev,
                       ppd.conv_crit.get() ? ppd.conv_crit->getVectorNormType()
                                           : MathLib::VecNormType::NORM2)