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       15454 : PiecewiseMultilinear::validParams()
      17             : {
      18       15454 :   InputParameters params = PiecewiseMultiInterpolation::validParams();
      19       30908 :   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       30908 :   params.addParam<Real>(
      25       30908 :       "epsilon", 1e-12, "Finite differencing parameter for gradient and time derivative");
      26       15454 :   return params;
      27           0 : }
      28             : 
      29         259 : PiecewiseMultilinear::PiecewiseMultilinear(const InputParameters & parameters)
      30         483 :   : PiecewiseMultiInterpolation(parameters), _epsilon(getParam<Real>("epsilon"))
      31             : {
      32         224 : }
      33             : 
      34             : Real
      35     1358972 : PiecewiseMultilinear::sample(const GridPoint & pt) const
      36             : {
      37     1358972 :   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     1358972 : PiecewiseMultilinear::sampleInternal(const Moose::GenericType<GridPoint, is_ad> pt) const
      49             : {
      50             :   using std::abs;
      51             :   /*
      52             :    * left contains the indices of the point to the 'left', 'down', etc, of pt
      53             :    * right contains the indices of the point to the 'right', 'up', etc, of pt
      54             :    * Hence, left and right define the vertices of the hypercube containing pt
      55             :    */
      56     1358972 :   GridIndex left(_dim);
      57     1358972 :   GridIndex right(_dim);
      58     5427670 :   for (unsigned int i = 0; i < _dim; ++i)
      59     4068698 :     getNeighborIndices(_grid[i], MetaPhysicL::raw_value(pt[i]), left[i], right[i]);
      60             : 
      61             :   /*
      62             :    * The following just loops through all the vertices of the
      63             :    * hypercube containing pt, evaluating the function at all
      64             :    * those vertices, and weighting the contributions to the
      65             :    * final result depending on the distance of pt from the vertex
      66             :    */
      67     1358972 :   Moose::GenericType<Real, is_ad> f = 0;
      68           0 :   Moose::GenericType<Real, is_ad> weight;
      69     1358972 :   GridIndex arg(_dim);
      70             :   // number of points in hypercube = 2^_dim
      71    12229684 :   for (unsigned int i = 0; i < (1u << _dim); ++i)
      72             :   {
      73    10870712 :     weight = 1;
      74    43533248 :     for (unsigned int j = 0; j < _dim; ++j)
      75             :       // shift i j-bits to the right and see if the result has a 0 as its right-most bit
      76    32662536 :       if ((i >> j) % 2 == 0)
      77             :       {
      78    16331268 :         arg[j] = left[j];
      79    16331268 :         if (left[j] != right[j])
      80     8698608 :           weight *= abs(pt[j] - _grid[j][right[j]]);
      81             :         else
      82             :           // unusual "end condition" case. weight by 0.5 because we will encounter this twice
      83     7632660 :           weight *= 0.5;
      84             :       }
      85             :       else
      86             :       {
      87    16331268 :         arg[j] = right[j];
      88    16331268 :         if (left[j] != right[j])
      89     8698608 :           weight *= abs(pt[j] - _grid[j][left[j]]);
      90             :         else
      91             :           // unusual "end condition" case. weight by 0.5 because we will encounter this twice
      92     7632660 :           weight *= 0.5;
      93             :       }
      94    10870712 :     f += _gridded_data->evaluateFcn(arg) * weight;
      95             :   }
      96             : 
      97             :   /*
      98             :    * finally divide by the volume of the hypercube
      99             :    */
     100     1358972 :   weight = 1;
     101     5427670 :   for (unsigned int dim = 0; dim < pt.size(); ++dim)
     102     4068698 :     if (left[dim] != right[dim])
     103     2170678 :       weight *= _grid[dim][right[dim]] - _grid[dim][left[dim]];
     104             :     else
     105             :       // unusual "end condition" case. weight by 1 to cancel the two 0.5 encountered previously
     106     1898020 :       weight *= 1;
     107             : 
     108     1358972 :   return f / weight;
     109           0 : }
     110             : 
     111             : RealGradient
     112           0 : PiecewiseMultilinear::gradient(Real t, const Point & p) const
     113             : {
     114           0 :   RealGradient grad;
     115             : 
     116             :   // sample center point
     117           0 :   auto s1 = sample(pointInGrid<false>(t, p));
     118             : 
     119             :   // sample epsilon steps in all directions
     120           0 :   for (const auto dir : _axes)
     121           0 :     if (dir < 3)
     122             :     {
     123           0 :       Point pp = p;
     124           0 :       pp(dir) += _epsilon;
     125           0 :       grad(dir) = (sample(pointInGrid<false>(t, pp)) - s1) / _epsilon;
     126             :     }
     127             : 
     128           0 :   return grad;
     129             : }
     130             : 
     131             : Real
     132           0 : PiecewiseMultilinear::timeDerivative(Real t, const Point & p) const
     133             : {
     134           0 :   return (sample(pointInGrid<false>(t + _epsilon, p)) - sample(pointInGrid<false>(t, p))) /
     135           0 :          _epsilon;
     136             : }