Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "LinearInterpolation.h"
      11             : #include "MooseUtils.h"
      12             : 
      13             : #include "ChainedReal.h"
      14             : 
      15             : #include <cassert>
      16             : #include <fstream>
      17             : #include <stdexcept>
      18             : 
      19        3422 : LinearInterpolation::LinearInterpolation(const std::vector<Real> & x,
      20             :                                          const std::vector<Real> & y,
      21        3422 :                                          const bool extrap)
      22        3422 :   : _x(x), _y(y), _extrap(extrap)
      23             : {
      24        3422 :   errorCheck();
      25        3434 : }
      26             : 
      27             : void
      28        3522 : LinearInterpolation::errorCheck()
      29             : {
      30        3522 :   if (_x.size() != _y.size())
      31           8 :     throw std::domain_error("Vectors are not the same length");
      32             : 
      33      121302 :   for (unsigned int i = 0; i + 1 < _x.size(); ++i)
      34      117792 :     if (_x[i] >= _x[i + 1])
      35             :     {
      36           4 :       std::ostringstream oss;
      37           4 :       oss << "x-values are not strictly increasing: x[" << i << "]: " << _x[i] << " x[" << i + 1
      38           4 :           << "]: " << _x[i + 1];
      39           4 :       throw std::domain_error(oss.str());
      40           4 :     }
      41        3510 : }
      42             : 
      43             : template <typename T>
      44             : T
      45      671419 : LinearInterpolation::sample(const T & x) const
      46             : {
      47             :   // this is a hard error
      48      671419 :   if (_x.empty())
      49           0 :     mooseError("Trying to evaluate an empty LinearInterpolation");
      50             : 
      51             :   // special case for single function nodes
      52      671419 :   if (_x.size() == 1)
      53           0 :     return _y[0];
      54             : 
      55             :   // endpoint cases
      56      671419 :   if (_extrap)
      57             :   {
      58        1200 :     if (x < _x[0])
      59         776 :       return _y[0] + (x - _x[0]) / (_x[1] - _x[0]) * (_y[1] - _y[0]);
      60         424 :     if (x >= _x.back())
      61         208 :       return _y.back() +
      62         304 :              (x - _x.back()) / (_x[_x.size() - 2] - _x.back()) * (_y[_y.size() - 2] - _y.back());
      63             :   }
      64             :   else
      65             :   {
      66      670219 :     if (x < _x[0])
      67         844 :       return _y[0];
      68      669375 :     if (x >= _x.back())
      69       20150 :       return _y.back();
      70             :   }
      71             : 
      72      649441 :   auto upper = std::upper_bound(_x.begin(), _x.end(), x);
      73      649441 :   const auto i = cast_int<std::size_t>(std::distance(_x.begin(), upper) - 1);
      74      649441 :   if (i == cast_int<std::size_t>(_x.size() - 1))
      75             :     // std::upper_bound returns the end() iterator if there are no elements that are
      76             :     // an upper bound to the value. Since x >= _x.back() has already returned above,
      77             :     // this means x is a NaN, so we return a NaN here.
      78           1 :     return std::nan("");
      79             :   else
      80      649440 :     return _y[i] + (_y[i + 1] - _y[i]) * (x - _x[i]) / (_x[i + 1] - _x[i]);
      81             : }
      82             : 
      83             : template Real LinearInterpolation::sample<Real>(const Real &) const;
      84             : template ADReal LinearInterpolation::sample<ADReal>(const ADReal &) const;
      85             : template ChainedReal LinearInterpolation::sample<ChainedReal>(const ChainedReal &) const;
      86             : 
      87             : template <typename T>
      88             : T
      89        4357 : LinearInterpolation::sampleDerivative(const T & x) const
      90             : {
      91             :   // endpoint cases
      92        4357 :   if (_extrap)
      93             :   {
      94         576 :     if (x <= _x[0])
      95         384 :       return (_y[1] - _y[0]) / (_x[1] - _x[0]);
      96         192 :     if (x >= _x.back())
      97          96 :       return (_y[_y.size() - 2] - _y.back()) / (_x[_x.size() - 2] - _x.back());
      98             :   }
      99             :   else
     100             :   {
     101        3781 :     if (x < _x[0])
     102           1 :       return 0.0;
     103        3780 :     if (x >= _x.back())
     104           0 :       return 0.0;
     105             :   }
     106             : 
     107        3876 :   auto upper = std::upper_bound(_x.begin(), _x.end(), x);
     108        3876 :   const auto i = cast_int<std::size_t>(std::distance(_x.begin(), upper) - 1);
     109        3876 :   if (i == cast_int<std::size_t>(_x.size() - 1))
     110             :     // std::upper_bound returns the end() iterator if there are no elements that are
     111             :     // an upper bound to the value. Since x >= _x.back() has already returned above,
     112             :     // this means x is a NaN, so we return a NaN here.
     113           1 :     return std::nan("");
     114             :   else
     115        3875 :     return (_y[i + 1] - _y[i]) / (_x[i + 1] - _x[i]);
     116             : }
     117             : 
     118             : template Real LinearInterpolation::sampleDerivative<Real>(const Real &) const;
     119             : template ADReal LinearInterpolation::sampleDerivative<ADReal>(const ADReal &) const;
     120             : template ChainedReal LinearInterpolation::sampleDerivative<ChainedReal>(const ChainedReal &) const;
     121             : 
     122             : Real
     123           0 : LinearInterpolation::integrate()
     124             : {
     125           0 :   Real answer = 0;
     126           0 :   for (unsigned int i = 1; i < _x.size(); ++i)
     127           0 :     answer += 0.5 * (_y[i] + _y[i - 1]) * (_x[i] - _x[i - 1]);
     128             : 
     129           0 :   return answer;
     130             : }
     131             : 
     132             : Real
     133           5 : LinearInterpolation::integratePartial(Real xA, Real xB) const
     134             : {
     135             :   // integral computation below will assume that x2 > x1; if this is not the
     136             :   // case, compute as if it is and then use identity to convert
     137             :   Real x1, x2;
     138             :   bool switch_bounds;
     139           5 :   if (MooseUtils::absoluteFuzzyEqual(xA, xB))
     140           0 :     return 0.0;
     141           5 :   else if (xB > xA)
     142             :   {
     143           4 :     x1 = xA;
     144           4 :     x2 = xB;
     145           4 :     switch_bounds = false;
     146             :   }
     147             :   else
     148             :   {
     149           1 :     x1 = xB;
     150           1 :     x2 = xA;
     151           1 :     switch_bounds = true;
     152             :   }
     153             : 
     154             :   // compute integral with knowledge that x2 > x1
     155           5 :   Real integral = 0.0;
     156             :   // find minimum i : x[i] > x; if x > x[n-1], i = n
     157           5 :   auto n = _x.size();
     158             :   const unsigned int i1 =
     159           5 :       x1 <= _x[n - 1] ? std::distance(_x.begin(), std::upper_bound(_x.begin(), _x.end(), x1)) : n;
     160             :   const unsigned int i2 =
     161           5 :       x2 <= _x[n - 1] ? std::distance(_x.begin(), std::upper_bound(_x.begin(), _x.end(), x2)) : n;
     162           5 :   unsigned int i = i1;
     163          22 :   while (i <= i2)
     164             :   {
     165          17 :     if (i == 0)
     166             :     {
     167             :       // note i1 = i
     168             :       Real integral1, integral2;
     169           2 :       if (_extrap)
     170             :       {
     171           0 :         const Real dydx = (_y[1] - _y[0]) / (_x[1] - _x[0]);
     172           0 :         const Real y1 = _y[0] + dydx * (x1 - _x[0]);
     173           0 :         integral1 = 0.5 * (y1 + _y[0]) * (_x[0] - x1);
     174           0 :         if (i2 == i)
     175             :         {
     176           0 :           const Real y2 = _y[0] + dydx * (x2 - _x[0]);
     177           0 :           integral2 = 0.5 * (y2 + _y[0]) * (_x[0] - x2);
     178             :         }
     179             :         else
     180           0 :           integral2 = 0.0;
     181             :       }
     182             :       else
     183             :       {
     184           2 :         integral1 = _y[0] * (_x[0] - x1);
     185           2 :         if (i2 == i)
     186           0 :           integral2 = _y[0] * (_x[0] - x2);
     187             :         else
     188           2 :           integral2 = 0.0;
     189             :       }
     190             : 
     191           2 :       integral += integral1 - integral2;
     192             :     }
     193          15 :     else if (i == n)
     194             :     {
     195             :       // note i2 = i
     196             :       Real integral1, integral2;
     197           2 :       if (_extrap)
     198             :       {
     199           0 :         const Real dydx = (_y[n - 1] - _y[n - 2]) / (_x[n - 1] - _x[n - 2]);
     200           0 :         const Real y2 = _y[n - 1] + dydx * (x2 - _x[n - 1]);
     201           0 :         integral2 = 0.5 * (y2 + _y[n - 1]) * (x2 - _x[n - 1]);
     202           0 :         if (i1 == n)
     203             :         {
     204           0 :           const Real y1 = _y[n - 1] + dydx * (x1 - _x[n - 1]);
     205           0 :           integral1 = 0.5 * (y1 + _y[n - 1]) * (x1 - _x[n - 1]);
     206             :         }
     207             :         else
     208           0 :           integral1 = 0.0;
     209             :       }
     210             :       else
     211             :       {
     212           2 :         integral2 = _y[n - 1] * (x2 - _x[n - 1]);
     213           2 :         if (i1 == n)
     214           0 :           integral1 = _y[n - 1] * (x1 - _x[n - 1]);
     215             :         else
     216           2 :           integral1 = 0.0;
     217             :       }
     218             : 
     219           2 :       integral += integral2 - integral1;
     220             :     }
     221             :     else
     222             :     {
     223             :       Real integral1;
     224          13 :       if (i == i1)
     225             :       {
     226           3 :         const Real dydx = (_y[i] - _y[i - 1]) / (_x[i] - _x[i - 1]);
     227           3 :         const Real y1 = _y[i - 1] + dydx * (x1 - _x[i - 1]);
     228           3 :         integral1 = 0.5 * (y1 + _y[i - 1]) * (x1 - _x[i - 1]);
     229             :       }
     230             :       else
     231          10 :         integral1 = 0.0;
     232             : 
     233             :       Real integral2;
     234          13 :       if (i == i2)
     235             :       {
     236           3 :         const Real dydx = (_y[i] - _y[i - 1]) / (_x[i] - _x[i - 1]);
     237           3 :         const Real y2 = _y[i - 1] + dydx * (x2 - _x[i - 1]);
     238           3 :         integral2 = 0.5 * (y2 + _y[i - 1]) * (x2 - _x[i - 1]);
     239             :       }
     240             :       else
     241          10 :         integral2 = 0.5 * (_y[i] + _y[i - 1]) * (_x[i] - _x[i - 1]);
     242             : 
     243          13 :       integral += integral2 - integral1;
     244             :     }
     245             : 
     246          17 :     i++;
     247             :   }
     248             : 
     249             :   // apply identity if bounds were switched
     250           5 :   if (switch_bounds)
     251           1 :     return -1.0 * integral;
     252             :   else
     253           4 :     return integral;
     254             : }
     255             : 
     256             : Real
     257           0 : LinearInterpolation::domain(int i) const
     258             : {
     259           0 :   return _x[i];
     260             : }
     261             : 
     262             : Real
     263           0 : LinearInterpolation::range(int i) const
     264             : {
     265           0 :   return _y[i];
     266             : }
     267             : 
     268             : unsigned int
     269         395 : LinearInterpolation::getSampleSize() const
     270             : {
     271         395 :   return _x.size();
     272             : }