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 "PiecewiseMultilinear.h"
      11             : #include "GriddedData.h"
      12             : 
      13             : registerMooseObject("MooseApp", PiecewiseMultilinear);
      14             : 
      15             : InputParameters
      16       14700 : PiecewiseMultilinear::validParams()
      17             : {
      18       14700 :   InputParameters params = PiecewiseMultiInterpolation::validParams();
      19       14700 :   params.addClassDescription(
      20             :       "PiecewiseMultilinear performs linear interpolation on 1D, 2D, 3D or 4D "
      21             :       "data.  The data_file specifies the axes directions and the function "
      22             :       "values.  If a point lies outside the data range, the appropriate end "
      23             :       "value is used.");
      24       44100 :   params.addParam<Real>(
      25       29400 :       "epsilon", 1e-12, "Finite differencing parameter for gradient and time derivative");
      26       14700 :   return params;
      27           0 : }
      28             : 
      29         243 : PiecewiseMultilinear::PiecewiseMultilinear(const InputParameters & parameters)
      30         243 :   : PiecewiseMultiInterpolation(parameters), _epsilon(getParam<Real>("epsilon"))
      31             : {
      32         208 : }
      33             : 
      34             : Real
      35     1200976 : PiecewiseMultilinear::sample(const GridPoint & pt) const
      36             : {
      37     1200976 :   return sampleInternal<false>(pt);
      38             : }
      39             : 
      40             : ADReal
      41           0 : PiecewiseMultilinear::sample(const ADGridPoint & pt) const
      42             : {
      43           0 :   return sampleInternal<true>(pt);
      44             : }
      45             : 
      46             : template <bool is_ad>
      47             : Moose::GenericType<Real, is_ad>
      48     1200976 : PiecewiseMultilinear::sampleInternal(const Moose::GenericType<GridPoint, is_ad> pt) const
      49             : {
      50             :   /*
      51             :    * left contains the indices of the point to the 'left', 'down', etc, of pt
      52             :    * right contains the indices of the point to the 'right', 'up', etc, of pt
      53             :    * Hence, left and right define the vertices of the hypercube containing pt
      54             :    */
      55     1200976 :   GridIndex left(_dim);
      56     1200976 :   GridIndex right(_dim);
      57     4797136 :   for (unsigned int i = 0; i < _dim; ++i)
      58     3596160 :     getNeighborIndices(_grid[i], MetaPhysicL::raw_value(pt[i]), left[i], right[i]);
      59             : 
      60             :   /*
      61             :    * The following just loops through all the vertices of the
      62             :    * hypercube containing pt, evaluating the function at all
      63             :    * those vertices, and weighting the contributions to the
      64             :    * final result depending on the distance of pt from the vertex
      65             :    */
      66     1200976 :   Moose::GenericType<Real, is_ad> f = 0;
      67           0 :   Moose::GenericType<Real, is_ad> weight;
      68     1200976 :   GridIndex arg(_dim);
      69             :   // number of points in hypercube = 2^_dim
      70    10809296 :   for (unsigned int i = 0; i < (1u << _dim); ++i)
      71             :   {
      72     9608320 :     weight = 1;
      73    38478656 :     for (unsigned int j = 0; j < _dim; ++j)
      74             :       // shift i j-bits to the right and see if the result has a 0 as its right-most bit
      75    28870336 :       if ((i >> j) % 2 == 0)
      76             :       {
      77    14435168 :         arg[j] = left[j];
      78    14435168 :         if (left[j] != right[j])
      79     7681920 :           weight *= std::abs(pt[j] - _grid[j][right[j]]);
      80             :         else
      81             :           // unusual "end condition" case. weight by 0.5 because we will encounter this twice
      82     6753248 :           weight *= 0.5;
      83             :       }
      84             :       else
      85             :       {
      86    14435168 :         arg[j] = right[j];
      87    14435168 :         if (left[j] != right[j])
      88     7681920 :           weight *= std::abs(pt[j] - _grid[j][left[j]]);
      89             :         else
      90             :           // unusual "end condition" case. weight by 0.5 because we will encounter this twice
      91     6753248 :           weight *= 0.5;
      92             :       }
      93     9608320 :     f += _gridded_data->evaluateFcn(arg) * weight;
      94             :   }
      95             : 
      96             :   /*
      97             :    * finally divide by the volume of the hypercube
      98             :    */
      99     1200976 :   weight = 1;
     100     4797136 :   for (unsigned int dim = 0; dim < pt.size(); ++dim)
     101     3596160 :     if (left[dim] != right[dim])
     102     1916912 :       weight *= _grid[dim][right[dim]] - _grid[dim][left[dim]];
     103             :     else
     104             :       // unusual "end condition" case. weight by 1 to cancel the two 0.5 encountered previously
     105     1679248 :       weight *= 1;
     106             : 
     107     1200976 :   return f / weight;
     108           0 : }
     109             : 
     110             : RealGradient
     111           0 : PiecewiseMultilinear::gradient(Real t, const Point & p) const
     112             : {
     113           0 :   RealGradient grad;
     114             : 
     115             :   // sample center point
     116           0 :   auto s1 = sample(pointInGrid<false>(t, p));
     117             : 
     118             :   // sample epsilon steps in all directions
     119           0 :   for (const auto dir : _axes)
     120           0 :     if (dir < 3)
     121             :     {
     122           0 :       Point pp = p;
     123           0 :       pp(dir) += _epsilon;
     124           0 :       grad(dir) = (sample(pointInGrid<false>(t, pp)) - s1) / _epsilon;
     125             :     }
     126             : 
     127           0 :   return grad;
     128             : }
     129             : 
     130             : Real
     131           0 : PiecewiseMultilinear::timeDerivative(Real t, const Point & p) const
     132             : {
     133           0 :   return (sample(pointInGrid<false>(t + _epsilon, p)) - sample(pointInGrid<false>(t, p))) /
     134           0 :          _epsilon;
     135             : }