Merge branch 'MoveCodeFromNumLib2MathLib' into 'master'

Code cleanup: Move code from NumLib to MathLib

See merge request ogs/ogs!4833

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

NumLib/ODESolver/TimeDiscretization.cpp

@@ -15,39 +15,6 @@
 
 namespace NumLib
 {
-double computeRelativeNorm(GlobalVector const& x,
-                           GlobalVector const& x_prev,
-                           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 - x_prev --> x - dx --> dx
-    MathLib::LinAlg::axpy(dx, -1.0, x_prev);
-    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
 {

NumLib/ODESolver/TimeDiscretization.h

@@ -18,19 +18,6 @@
 
 namespace NumLib
 {
-/**
- * Compute and return the relative norm of solutions between two
- * successive time steps by \f$ \|u^{n+1}-u^{n}\|/\|u^{n+1}\| \f$.
- *
- * \param x         The current solution
- * \param x_prev    The previous solution
- * \param norm_type The norm type of global vector
- * \return          \f$ \|u^{n+1}-u^{n}\|/\|u^{n+1}\| \f$.
- */
-double computeRelativeNorm(GlobalVector const& x,
-                           GlobalVector const& x_prev,
-                           MathLib::VecNormType norm_type);
-
 //! \addtogroup ODESolver
 //! @{
 

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)