From 7c92de16a57dde5c26056aa5a1efa3a5e931a443 Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Thu, 7 Dec 2023 14:13:56 +0100
Subject: [PATCH] [NL|MathLib] Mv computeRelativeNorm from NumLib to MathLib
---
MathLib/LinAlg/LinAlg.h | 45 +++++++++++++++++++++++++
NumLib/ODESolver/TimeDiscretization.cpp | 34 -------------------
NumLib/ODESolver/TimeDiscretization.h | 11 ------
ProcessLib/TimeLoop.cpp | 2 +-
4 files changed, 46 insertions(+), 46 deletions(-)
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)
--
GitLab