/**
 * \copyright
 * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
 *            Distributed under a Modified BSD License.
 *              See accompanying file LICENSE.txt or
 *              http://www.opengeosys.org/project/license
 *
 *  \file
 *  Created on March 31, 2017, 10:30 AM
 */

#include "TimeDiscretization.h"

#include "MathLib/LinAlg/MatrixVectorTraits.h"

namespace NumLib
{
double TimeDiscretization::computeRelativeChangeFromPreviousTimestep(
    GlobalVector const& x,
    GlobalVector const& x_old,
    MathLib::VecNormType norm_type)
{
    if (norm_type == MathLib::VecNormType::INVALID)
    {
        OGS_FATAL("An invalid norm type has been passed");
    }

    if (!_dx)
    {
        _dx = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(x);
    }

    auto& dx = *_dx;
    MathLib::LinAlg::copy(x, dx);  // copy x to dx.

    // dx = x - x_old --> x - dx --> dx
    MathLib::LinAlg::axpy(dx, -1.0, x_old);
    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();
}
}  // end of namespace NumLib