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        3496 : LinearInterpolation::LinearInterpolation(const std::vector<Real> & x,
      20             :                                          const std::vector<Real> & y,
      21        3496 :                                          const bool extrap)
      22        3496 :   : _x(x), _y(y), _extrap(extrap)
      23             : {
      24        3496 :   errorCheck();
      25        3505 : }
      26             : 
      27             : void
      28        3596 : LinearInterpolation::errorCheck()
      29             : {
      30        3596 :   if (_x.size() != _y.size())
      31           6 :     throw std::domain_error("Vectors are not the same length");
      32             : 
      33      141225 :   for (unsigned int i = 0; i + 1 < _x.size(); ++i)
      34      137638 :     if (_x[i] >= _x[i + 1])
      35             :     {
      36           3 :       std::ostringstream oss;
      37           3 :       oss << "x-values are not strictly increasing: x[" << i << "]: " << _x[i] << " x[" << i + 1
      38           3 :           << "]: " << _x[i + 1];
      39           3 :       throw std::domain_error(oss.str());
      40           3 :     }
      41        3587 : }
      42             : 
      43             : template <typename T>
      44             : T
      45      679749 : LinearInterpolation::sample(const T & x) const
      46             : {
      47             :   // this is a hard error
      48      679749 :   if (_x.empty())
      49           0 :     mooseError("Trying to evaluate an empty LinearInterpolation");
      50             : 
      51             :   // special case for single function nodes
      52      679749 :   if (_x.size() == 1)
      53           0 :     return _y[0];
      54             : 
      55             :   // endpoint cases
      56      679749 :   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      678549 :     if (x < _x[0])
      67        1390 :       return _y[0];
      68      677159 :     if (x >= _x.back())
      69       20378 :       return _y.back();
      70             :   }
      71             : 
      72      656997 :   auto upper = std::upper_bound(_x.begin(), _x.end(), x);
      73     1313994 :   const auto i = cast_int<std::size_t>(std::distance(_x.begin(), upper) - 1);
      74      656997 :   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           2 :     return std::nan("");
      79             :   else
      80      656995 :     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        4386 : LinearInterpolation::sampleDerivative(const T & x) const
      90             : {
      91             :   // endpoint cases
      92        4386 :   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        3810 :     if (x < _x[0])
     102           2 :       return 0.0;
     103        3808 :     if (x >= _x.back())
     104           0 :       return 0.0;
     105             :   }
     106             : 
     107        3904 :   auto upper = std::upper_bound(_x.begin(), _x.end(), x);
     108        7808 :   const auto i = cast_int<std::size_t>(std::distance(_x.begin(), upper) - 1);
     109        3904 :   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           2 :     return std::nan("");
     114             :   else
     115        3902 :     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             : template <typename T>
     133             : T
     134          10 : LinearInterpolation::integratePartial(const T & xA, const T & xB) const
     135             : {
     136             :   // integral computation below will assume that x2 > x1; if this is not the
     137             :   // case, compute as if it is and then use identity to convert
     138             :   const T *x1_ptr, *x2_ptr;
     139             :   bool switch_bounds;
     140          10 :   const Real rawA = MetaPhysicL::raw_value(xA), rawB = MetaPhysicL::raw_value(xB);
     141             :   Real raw1, raw2;
     142          10 :   if (MooseUtils::absoluteFuzzyEqual(rawA, rawB))
     143           0 :     return 0.0;
     144          10 :   else if (rawB > rawA)
     145             :   {
     146           8 :     x1_ptr = &xA;
     147           8 :     x2_ptr = &xB;
     148           8 :     raw1 = rawA;
     149           8 :     raw2 = rawB;
     150           8 :     switch_bounds = false;
     151             :   }
     152             :   else
     153             :   {
     154           2 :     x1_ptr = &xB;
     155           2 :     x2_ptr = &xA;
     156           2 :     raw1 = rawB;
     157           2 :     raw2 = rawA;
     158           2 :     switch_bounds = true;
     159             :   }
     160          10 :   const auto & x1 = *x1_ptr;
     161          10 :   const auto & x2 = *x2_ptr;
     162             : 
     163             :   // compute integral with knowledge that x2 > x1
     164          10 :   T integral = 0.0;
     165             :   // find minimum i : x[i] > x; if x > x[n-1], i = n
     166          10 :   const auto n = _x.size();
     167           0 :   const unsigned int i1 =
     168          20 :       raw1 <= _x[n - 1] ? std::distance(_x.begin(), std::upper_bound(_x.begin(), _x.end(), raw1))
     169             :                         : n;
     170           4 :   const unsigned int i2 =
     171          16 :       raw2 <= _x[n - 1] ? std::distance(_x.begin(), std::upper_bound(_x.begin(), _x.end(), raw2))
     172             :                         : n;
     173          10 :   unsigned int i = i1;
     174          44 :   while (i <= i2)
     175             :   {
     176          34 :     if (i == 0)
     177             :     {
     178             :       // note i1 = i
     179           0 :       T integral1, integral2;
     180           4 :       if (_extrap)
     181             :       {
     182           0 :         const auto dydx = (_y[1] - _y[0]) / (_x[1] - _x[0]);
     183           0 :         const auto y1 = _y[0] + dydx * (x1 - _x[0]);
     184           0 :         integral1 = 0.5 * (y1 + _y[0]) * (_x[0] - x1);
     185           0 :         if (i2 == i)
     186             :         {
     187           0 :           const auto y2 = _y[0] + dydx * (x2 - _x[0]);
     188           0 :           integral2 = 0.5 * (y2 + _y[0]) * (_x[0] - x2);
     189           0 :         }
     190             :         else
     191           0 :           integral2 = 0.0;
     192           0 :       }
     193             :       else
     194             :       {
     195           4 :         integral1 = _y[0] * (_x[0] - x1);
     196           4 :         if (i2 == i)
     197           0 :           integral2 = _y[0] * (_x[0] - x2);
     198             :         else
     199           4 :           integral2 = 0.0;
     200             :       }
     201             : 
     202           4 :       integral += integral1 - integral2;
     203           0 :     }
     204          30 :     else if (i == n)
     205             :     {
     206             :       // note i2 = i
     207           0 :       T integral1, integral2;
     208           4 :       if (_extrap)
     209             :       {
     210           0 :         const auto dydx = (_y[n - 1] - _y[n - 2]) / (_x[n - 1] - _x[n - 2]);
     211           0 :         const auto y2 = _y[n - 1] + dydx * (x2 - _x[n - 1]);
     212           0 :         integral2 = 0.5 * (y2 + _y[n - 1]) * (x2 - _x[n - 1]);
     213           0 :         if (i1 == n)
     214             :         {
     215           0 :           const auto y1 = _y[n - 1] + dydx * (x1 - _x[n - 1]);
     216           0 :           integral1 = 0.5 * (y1 + _y[n - 1]) * (x1 - _x[n - 1]);
     217           0 :         }
     218             :         else
     219           0 :           integral1 = 0.0;
     220           0 :       }
     221             :       else
     222             :       {
     223           4 :         integral2 = _y[n - 1] * (x2 - _x[n - 1]);
     224           4 :         if (i1 == n)
     225           0 :           integral1 = _y[n - 1] * (x1 - _x[n - 1]);
     226             :         else
     227           4 :           integral1 = 0.0;
     228             :       }
     229             : 
     230           4 :       integral += integral2 - integral1;
     231           0 :     }
     232             :     else
     233             :     {
     234           0 :       T integral1;
     235          26 :       if (i == i1)
     236             :       {
     237           6 :         const auto dydx = (_y[i] - _y[i - 1]) / (_x[i] - _x[i - 1]);
     238           6 :         const auto y1 = _y[i - 1] + dydx * (x1 - _x[i - 1]);
     239           6 :         integral1 = 0.5 * (y1 + _y[i - 1]) * (x1 - _x[i - 1]);
     240           0 :       }
     241             :       else
     242          20 :         integral1 = 0.0;
     243             : 
     244           0 :       T integral2;
     245          26 :       if (i == i2)
     246             :       {
     247           6 :         const auto dydx = (_y[i] - _y[i - 1]) / (_x[i] - _x[i - 1]);
     248           6 :         const auto y2 = _y[i - 1] + dydx * (x2 - _x[i - 1]);
     249           6 :         integral2 = 0.5 * (y2 + _y[i - 1]) * (x2 - _x[i - 1]);
     250           0 :       }
     251             :       else
     252          20 :         integral2 = 0.5 * (_y[i] + _y[i - 1]) * (_x[i] - _x[i - 1]);
     253             : 
     254          26 :       integral += integral2 - integral1;
     255           0 :     }
     256             : 
     257          34 :     i++;
     258             :   }
     259             : 
     260             :   // apply identity if bounds were switched
     261          10 :   if (switch_bounds)
     262           2 :     return -1.0 * integral;
     263             :   else
     264           8 :     return integral;
     265           0 : }
     266             : 
     267             : template Real LinearInterpolation::integratePartial(const Real &, const Real &) const;
     268             : template ADReal LinearInterpolation::integratePartial(const ADReal &, const ADReal &) const;
     269             : 
     270             : Real
     271           0 : LinearInterpolation::domain(int i) const
     272             : {
     273           0 :   return _x[i];
     274             : }
     275             : 
     276             : Real
     277           0 : LinearInterpolation::range(int i) const
     278             : {
     279           0 :   return _y[i];
     280             : }
     281             : 
     282             : unsigned int
     283         402 : LinearInterpolation::getSampleSize() const
     284             : {
     285         402 :   return _x.size();
     286             : }