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 "MooseUtils.h"
      11             : #include "PiecewiseMultiInterpolation.h"
      12             : #include "GriddedData.h"
      13             : 
      14             : InputParameters
      15       43255 : PiecewiseMultiInterpolation::validParams()
      16             : {
      17       43255 :   InputParameters params = Function::validParams();
      18       43255 :   params.addParam<FileName>(
      19             :       "data_file",
      20             :       "File holding data for use with PiecewiseMultiInterpolation.  Format: any empty line and any "
      21             :       "line "
      22             :       "beginning with # are ignored, all other lines are assumed to contain relevant information.  "
      23             :       "The file must begin with specification of the grid.  This is done through lines containing "
      24             :       "the keywords: AXIS X; AXIS Y; AXIS Z; or AXIS T.  Immediately following the keyword line "
      25             :       "must be a space-separated line of real numbers which define the grid along the specified "
      26             :       "axis.  These data must be monotonically increasing.  After all the axes and their grids "
      27             :       "have been specified, there must be a line that is DATA.  Following that line, function "
      28             :       "values are given in the correct order (they may be on individual lines, or be "
      29             :       "space-separated on a number of lines).  When the function is evaluated, f[i,j,k,l] "
      30             :       "corresponds to the i + j*Ni + k*Ni*Nj + l*Ni*Nj*Nk data value.  Here i>=0 corresponding to "
      31             :       "the index along the first AXIS, j>=0 corresponding to the index along the second AXIS, etc, "
      32             :       "and Ni = number of grid points along the first AXIS, etc.");
      33       43255 :   return params;
      34           0 : }
      35             : 
      36         256 : PiecewiseMultiInterpolation::PiecewiseMultiInterpolation(const InputParameters & parameters)
      37             :   : Function(parameters),
      38         256 :     _gridded_data(std::make_unique<GriddedData>(getParam<FileName>("data_file"))),
      39         491 :     _dim(_gridded_data->getDim())
      40             : {
      41         235 :   _gridded_data->getAxes(_axes);
      42         235 :   _gridded_data->getGrid(_grid);
      43             : 
      44             :   // GriddedData does not require monotonicity of axes, but we do
      45         561 :   for (unsigned int i = 0; i < _dim; ++i)
      46        1621 :     for (unsigned int j = 1; j < _grid[i].size(); ++j)
      47        1295 :       if (_grid[i][j - 1] >= _grid[i][j])
      48           7 :         mooseError("PiecewiseMultiInterpolation needs monotonically-increasing axis data.  Axis ",
      49             :                    i,
      50             :                    " contains non-monotonicity at value ",
      51           7 :                    _grid[i][j]);
      52             : 
      53             :   // GriddedData does not demand that each axis is independent, but we do
      54         228 :   std::set<int> s(_axes.begin(), _axes.end());
      55         228 :   if (s.size() != _dim)
      56           7 :     mooseError(
      57             :         "PiecewiseMultiInterpolation needs the AXES to be independent.  Check the AXIS lines in "
      58             :         "your data file.");
      59         221 : }
      60             : 
      61         221 : PiecewiseMultiInterpolation::~PiecewiseMultiInterpolation() {}
      62             : 
      63             : template <bool is_ad>
      64             : Moose::GenericType<PiecewiseMultiInterpolation::GridPoint, is_ad>
      65     1202000 : PiecewiseMultiInterpolation::pointInGrid(const Moose::GenericType<Real, is_ad> & t,
      66             :                                          const Moose::GenericType<Point, is_ad> & p) const
      67             : {
      68             :   // convert the inputs to an input to the sample function using _axes
      69     1202000 :   Moose::GenericType<GridPoint, is_ad> point_in_grid(_dim);
      70     4800208 :   for (unsigned int i = 0; i < _dim; ++i)
      71             :   {
      72     3598208 :     if (_axes[i] < 3)
      73     2401184 :       point_in_grid[i] = p(_axes[i]);
      74     1197024 :     else if (_axes[i] == 3) // the time direction
      75     1197024 :       point_in_grid[i] = t;
      76             :   }
      77     1202000 :   return point_in_grid;
      78           0 : }
      79             : 
      80             : template PiecewiseMultiInterpolation::GridPoint
      81             : PiecewiseMultiInterpolation::pointInGrid<false>(const Real &, const Point &) const;
      82             : template PiecewiseMultiInterpolation::ADGridPoint
      83             : PiecewiseMultiInterpolation::pointInGrid<true>(const ADReal &, const ADPoint &) const;
      84             : 
      85             : Real
      86     1202000 : PiecewiseMultiInterpolation::value(Real t, const Point & p) const
      87             : {
      88     1202000 :   return sample(pointInGrid<false>(t, p));
      89             : }
      90             : 
      91             : ADReal
      92           0 : PiecewiseMultiInterpolation::value(const ADReal & t, const ADPoint & p) const
      93             : {
      94           0 :   return sample(pointInGrid<true>(t, p));
      95             : }
      96             : 
      97             : ADReal
      98           0 : PiecewiseMultiInterpolation::sample(const ADGridPoint &) const
      99             : {
     100           0 :   mooseError("The AD variant of 'sample' needs to be implemented");
     101             : }
     102             : 
     103             : void
     104     3598208 : PiecewiseMultiInterpolation::getNeighborIndices(std::vector<Real> in_arr,
     105             :                                                 Real x,
     106             :                                                 unsigned int & lower_x,
     107             :                                                 unsigned int & upper_x) const
     108             : {
     109     3598208 :   int N = in_arr.size();
     110     3598208 :   if (x <= in_arr[0])
     111             :   {
     112      556352 :     lower_x = 0;
     113      556352 :     upper_x = 0;
     114             :   }
     115     3041856 :   else if (x >= in_arr[N - 1])
     116             :   {
     117      616272 :     lower_x = N - 1;
     118      616272 :     upper_x = N - 1;
     119             :   }
     120             :   else
     121             :   {
     122             :     // returns up which points at the first element in inArr that is not less than x
     123     2425584 :     std::vector<double>::iterator up = std::lower_bound(in_arr.begin(), in_arr.end(), x);
     124             : 
     125             :     // std::distance returns std::difference_type, which can be negative in theory, but
     126             :     // in this context will always be >=0.  Therefore the explicit cast is just to shut
     127             :     // the compiler up.
     128     2425584 :     upper_x = static_cast<unsigned int>(std::distance(in_arr.begin(), up));
     129     2425584 :     if (in_arr[upper_x] == x)
     130      507648 :       lower_x = upper_x;
     131             :     else
     132     1917936 :       lower_x = upper_x - 1;
     133             :   }
     134     3598208 : }