Skip to content

Commit

Permalink
[NL|MathLib] Mv computeRelativeNorm from NumLib to MathLib
Browse files Browse the repository at this point in the history
  • Loading branch information
TomFischer committed Dec 12, 2023
1 parent 8f51376 commit 7c92de1
Show file tree
Hide file tree
Showing 4 changed files with 46 additions and 46 deletions.
45 changes: 45 additions & 0 deletions MathLib/LinAlg/LinAlg.h
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
#pragma once

#include <cassert>

#include "BaseLib/Error.h"
#include "LinAlgEnums.h"

Expand Down Expand Up @@ -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
34 changes: 0 additions & 34 deletions NumLib/ODESolver/TimeDiscretization.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
{
Expand Down
11 changes: 0 additions & 11 deletions NumLib/ODESolver/TimeDiscretization.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
//! @{

Expand Down
2 changes: 1 addition & 1 deletion ProcessLib/TimeLoop.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down

0 comments on commit 7c92de1

Please sign in to comment.