Revision c565eaab258c7bb2fd33a91d03ae67d4f4edc75a authored by Christoph Lehmann on 19 September 2023, 05:43:08 UTC, committed by Christoph Lehmann on 19 September 2023, 05:43:08 UTC
Draft: RTP Acceleration Step 2: Compute linear solver only upon timestep change See merge request ogs/ogs!4730
PiecewiseLinearMonotonicCurve.cpp
/**
* \file
* \copyright
* Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*
* Created on November 11, 2016, 10:49 AM
*/
#include "PiecewiseLinearMonotonicCurve.h"
#include <algorithm>
#include <cmath>
#include <limits>
#include "BaseLib/Error.h"
namespace MathLib
{
PiecewiseLinearMonotonicCurve::PiecewiseLinearMonotonicCurve(
std::vector<double> x, std::vector<double> y)
: PiecewiseLinearInterpolation(std::move(x), std::move(y), false)
{
if (!isStrongMonotonic())
{
OGS_FATAL("The given curve is not strong monotonic.");
}
}
bool PiecewiseLinearMonotonicCurve::isStrongMonotonic() const
{
const double gradient0 = getDerivative(supp_pnts_[0]);
if (std::abs(gradient0) < std::numeric_limits<double>::min())
{
return false;
}
return std::none_of(supp_pnts_.begin(), supp_pnts_.end(),
[&](const double p)
{ return this->getDerivative(p) * gradient0 <= 0.; });
}
double PiecewiseLinearMonotonicCurve::getInverseVariable(const double y) const
{
std::size_t interval_idx = 0;
if (values_at_supp_pnts_.front() < values_at_supp_pnts_.back())
{
if (y <= values_at_supp_pnts_.front())
{
return supp_pnts_[0];
}
if (y >= values_at_supp_pnts_.back())
{
return supp_pnts_[supp_pnts_.size() - 1];
}
// search interval that has the point inside
auto const& it(std::lower_bound(values_at_supp_pnts_.begin(),
values_at_supp_pnts_.end(), y));
interval_idx = std::distance(values_at_supp_pnts_.begin(), it) - 1;
}
else
{
if (y >= values_at_supp_pnts_.front())
{
return supp_pnts_[0];
}
if (y <= values_at_supp_pnts_.back())
{
return supp_pnts_[supp_pnts_.size() - 1];
}
// search interval in the reverse direction for the point inside
auto const& it(std::lower_bound(values_at_supp_pnts_.rbegin(),
values_at_supp_pnts_.rend(), y));
interval_idx = values_at_supp_pnts_.size() -
(std::distance(values_at_supp_pnts_.rbegin(), it)) - 1;
}
const double xi_1 = supp_pnts_[interval_idx + 1];
const double xi = supp_pnts_[interval_idx];
const double yi_1 = values_at_supp_pnts_[interval_idx + 1];
const double yi = values_at_supp_pnts_[interval_idx];
// compute gradient: m = (x_{i+1} - x_i) / (y_{i+1} - y_i)
const double m = (xi_1 - xi) / (yi_1 - yi);
// compute the variable by linear interpolation: x = m * (y - y_i) + x_i,
// and then return the result.
return m * (y - yi) + xi;
}
} // namespace MathLib
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...