TableFunction.hpp

Go to the documentation of this file.

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2024 TotalEnergies
 * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2023-2024 Chevron
 * Copyright (c) 2019-     GEOS/GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */
 
#ifndef GEOS_FUNCTIONS_TABLEFUNCTION_HPP_
#define GEOS_FUNCTIONS_TABLEFUNCTION_HPP_
 
#include "FunctionBase.hpp"
 
#include "common/format/EnumStrings.hpp"
#include "LvArray/src/tensorOps.hpp"
#include "common/format/table/TableFormatter.hpp"
#include "common/Units.hpp"
 
namespace geos
{
 
class TableFunction : public FunctionBase
{
public:
 
  enum class InterpolationType : integer
  {
    Linear,
    Nearest,
    Upper,
    Lower
  };
 
  struct OutputOptions
  {
    bool writeCSV;
    bool writeInLog;
  };
 
  static constexpr integer maxDimensions = 4;
 
  class KernelWrapper
  {
public:
 
 
    KernelWrapper() = default;
    KernelWrapper( KernelWrapper const & ) = default;
    KernelWrapper( KernelWrapper && ) = default;
    KernelWrapper & operator=( KernelWrapper const & ) = default;
 
    KernelWrapper & operator=( KernelWrapper && other )
    {
      m_coordinates = std::move( other.m_coordinates );
      m_values = std::move( other.m_values );
      m_interpolationMethod = other.m_interpolationMethod;
      return *this;
    }
 
 
    template< typename IN_ARRAY >
    GEOS_HOST_DEVICE
    real64 compute( IN_ARRAY const & input ) const;
 
    template< typename IN_ARRAY, typename OUT_ARRAY >
    GEOS_HOST_DEVICE
    real64 compute( IN_ARRAY const & input, OUT_ARRAY && derivatives ) const;
 
    void move( LvArray::MemorySpace const space, bool const touch )
    {
      m_coordinates.move( space, touch );
      m_values.move( space, touch );
    }
 
private:
 
    friend class TableFunction; // Allow only parent class to construct the wrapper
 
    KernelWrapper( InterpolationType interpolationMethod,
                   ArrayOfArraysView< real64 const > const & coordinates,
                   arrayView1d< real64 const > const & values );
 
    template< typename IN_ARRAY >
    GEOS_HOST_DEVICE
    real64
    interpolateLinear( IN_ARRAY const & input ) const;
 
    template< typename IN_ARRAY, typename OUT_ARRAY >
    GEOS_HOST_DEVICE
    real64
    interpolateLinear( IN_ARRAY const & input, OUT_ARRAY && derivatives ) const;
 
    template< typename IN_ARRAY >
    GEOS_HOST_DEVICE
    real64
    interpolateRound( IN_ARRAY const & input ) const;
 
    template< typename IN_ARRAY >
    GEOS_HOST_DEVICE
    real64
    getCoord( IN_ARRAY const & input, localIndex dim, InterpolationType interpolationMethod ) const;
 
    template< typename IN_ARRAY, typename OUT_ARRAY >
    GEOS_HOST_DEVICE
    real64
    interpolateRound( IN_ARRAY const & input, OUT_ARRAY && derivatives ) const;
 
    TableFunction::InterpolationType m_interpolationMethod = InterpolationType::Linear;
 
    ArrayOfArraysView< real64 const > m_coordinates;
 
    arrayView1d< real64 const > m_values;
  };
 
  TableFunction( const string & name,
                 dataRepository::Group * const parent );
 
  static string catalogName() { return "TableFunction"; }
 
  virtual void initializeFunction() override;
 
  void reInitializeFunction();
 
  void initializePostSubGroups() override;
 
  virtual void evaluate( dataRepository::Group const & group,
                         real64 const time,
                         SortedArrayView< localIndex const > const & set,
                         arrayView1d< real64 > const & result ) const override final
  {
    FunctionBase::evaluateT< TableFunction, parallelHostPolicy >( group, time, set, result );
  }
 
  virtual real64 evaluate( real64 const * const input ) const override final;
 
  real64 getCoord( real64 const * const input, localIndex dim, InterpolationType interpolationMethod ) const;
 
  void checkCoord( real64 coord, localIndex dim ) const;
 
  integer numDimensions() const { return LvArray::integerConversion< integer >( m_coordinates.size() ); }
 
  ArrayOfArraysView< real64 const > getCoordinates() const { return m_coordinates.toViewConst(); }
 
  ArrayOfArraysView< real64 > getCoordinates() { return m_coordinates.toView(); }
 
  arrayView1d< real64 const > getValues() const { return m_values.toViewConst(); }
 
  array1d< real64 > & getValues() { return m_values; }
 
  InterpolationType getInterpolationMethod() const { return m_interpolationMethod; }
 
  units::Unit getDimUnit( localIndex const dim ) const
  {
    return size_t(dim) < m_dimUnits.size() ? m_dimUnits[dim] : units::Unknown;
  }
 
  void setInterpolationMethod( InterpolationType const method );
 
  void setTableCoordinates( array1d< real64_array > const & coordinates,
                            std::vector< units::Unit > const & dimUnits = {} );
 
  void setDimUnits( std::vector< units::Unit > const & dimUnits )
  {
    m_dimUnits = dimUnits;
  }
 
  void setTableValues( real64_array values, units::Unit unit = units::Unknown );
 
  void setValueUnits( units::Unit unit )
  {
    m_valueUnit = unit;
  }
 
  units::Unit getValueUnit() const { return m_valueUnit; }
 
  string getTableDescription() const;
 
  string getCoordsDescription( integer dimId, bool shortUnitsToVariables ) const;
 
  string getValuesDescription() const;
 
  void outputTableData( OutputOptions const outputOpts ) const;
 
  KernelWrapper createKernelWrapper() const;
 
  struct viewKeyStruct
  {
    static constexpr char const * coordinatesString() { return "coordinates"; }
    static constexpr char const * valuesString() { return "values"; }
    static constexpr char const * interpolationString() { return "interpolation"; }
    static constexpr char const * coordinateFilesString() { return "coordinateFiles"; }
    static constexpr char const * voxelFileString() { return "voxelFile"; }
    static constexpr char const * writeCSVFlagString() { return "writeCSV"; }
  };
 
private:
 
  void readFile( string const & filename, array1d< real64 > & target );
 
 
  array1d< real64 > m_tableCoordinates1D;
 
  path_array m_coordinateFiles;
 
  Path m_voxelFile;
 
  InterpolationType m_interpolationMethod;
 
  ArrayOfArrays< real64 > m_coordinates;
 
  array1d< real64 > m_values;
 
  std::vector< units::Unit > m_dimUnits;
 
  units::Unit m_valueUnit;
 
  KernelWrapper m_kernelWrapper;
 
  integer m_writeCSV;
};
template< typename IN_ARRAY >
GEOS_HOST_DEVICE
GEOS_FORCE_INLINE
real64
TableFunction::KernelWrapper::compute( IN_ARRAY const & input ) const
{
  if( m_interpolationMethod == TableFunction::InterpolationType::Linear )
  {
    return interpolateLinear( input );
  }
  else // Nearest, Upper, Lower interpolation methods
  {
    return interpolateRound( input );
  }
}
 
template< typename IN_ARRAY >
GEOS_HOST_DEVICE
GEOS_FORCE_INLINE
real64
TableFunction::KernelWrapper::interpolateLinear( IN_ARRAY const & input ) const
{
  integer const numDimensions = LvArray::integerConversion< integer >( m_coordinates.size() );
  localIndex bounds[maxDimensions][2]{};
  real64 weights[maxDimensions][2]{};
 
  // Determine position, weights
  for( localIndex dim = 0; dim < numDimensions; ++dim )
  {
    arraySlice1d< real64 const > const coords = m_coordinates[dim];
    if( input[dim] <= coords[0] )
    {
      // Coordinate is to the left of this axis
      bounds[dim][0] = 0;
      bounds[dim][1] = 0;
      weights[dim][0] = 0.0;
      weights[dim][1] = 1.0;
    }
    else if( input[dim] >= coords[coords.size() - 1] )
    {
      // Coordinate is to the right of this axis
      bounds[dim][0] = coords.size() - 1;
      bounds[dim][1] = bounds[dim][0];
      weights[dim][0] = 1.0;
      weights[dim][1] = 0.0;
    }
    else
    {
      // Find the coordinate index
      // Sergey's note: find uses a binary search...  If we assume coordinates are
      // evenly spaced, we can speed things up considerably
      // Mel's note: As we cannot be sure coords are evenly spaced,
      // - Either we insert coords to get even spacing ( /!\ memory consumption ),
      // - Or we can use an interpolation search with an hint array which would be linearly interpolated ( benchmark ).
      auto const lower = LvArray::sortedArrayManipulation::find( coords.begin(), coords.size(), input[dim] );
      bounds[dim][1] = LvArray::integerConversion< localIndex >( lower );
      bounds[dim][0] = bounds[dim][1] - 1;
 
      real64 const dx = coords[bounds[dim][1]] - coords[bounds[dim][0]];
      weights[dim][0] = 1.0 - ( input[dim] - coords[bounds[dim][0]]) / dx;
      weights[dim][1] = 1.0 - weights[dim][0];
    }
  }
 
  // Calculate the result
  real64 value = 0.0;
  integer const numCorners = 1 << numDimensions;
  for( integer point = 0; point < numCorners; ++point )
  {
    // Find array index
    localIndex tableIndex = 0;
    localIndex stride = 1;
    for( integer dim = 0; dim < numDimensions; ++dim )
    {
      integer const corner = (point >> dim) & 1;
      tableIndex += bounds[dim][corner] * stride;
      stride *= m_coordinates.sizeOfArray( dim );
    }
 
    // Determine weighted value
    real64 cornerValue = m_values[tableIndex];
    for( integer dim = 0; dim < numDimensions; ++dim )
    {
      integer const corner = (point >> dim) & 1;
      cornerValue *= weights[dim][corner];
    }
    value += cornerValue;
  }
  return value;
}
 
template< typename IN_ARRAY >
GEOS_HOST_DEVICE
GEOS_FORCE_INLINE
real64
TableFunction::KernelWrapper::interpolateRound( IN_ARRAY const & input ) const
{
  integer const numDimensions = LvArray::integerConversion< integer >( m_coordinates.size() );
 
  // Determine the index to the nearest table entry
  localIndex tableIndex = 0;
  localIndex stride = 1;
  for( integer dim = 0; dim < numDimensions; ++dim )
  {
    arraySlice1d< real64 const > const coords = m_coordinates[dim];
    // Determine the index along each table axis
    localIndex subIndex;
    if( input[dim] <= coords[0] )
    {
      // Coordinate is to the left of the table axis
      subIndex = 0;
    }
    else if( input[dim] >= coords[coords.size() - 1] )
    {
      // Coordinate is to the right of the table axis
      subIndex = coords.size() - 1;
    }
    else
    {
      // Coordinate is within the table axis
      // Note: find() will return the index of the upper table vertex
      auto const lower = LvArray::sortedArrayManipulation::find( coords.begin(), coords.size(), input[dim] );
      subIndex = LvArray::integerConversion< localIndex >( lower );
 
      // Interpolation types:
      //   - Nearest returns the value of the closest table vertex
      //   - Upper returns the value of the next table vertex
      //   - Lower returns the value of the previous table vertex
      if( m_interpolationMethod == TableFunction::InterpolationType::Nearest )
      {
        if( ( input[dim] - coords[subIndex - 1]) <= ( coords[subIndex] - input[dim]) )
        {
          --subIndex;
        }
      }
      else if( m_interpolationMethod == TableFunction::InterpolationType::Lower )
      {
        if( subIndex > 0 )
        {
          --subIndex;
        }
      }
    }
 
    // Increment the global table index
    tableIndex += subIndex * stride;
    stride *= coords.size();
  }
 
  // Retrieve the nearest value
  return m_values[tableIndex];
}
 
template< typename IN_ARRAY >
GEOS_HOST_DEVICE
GEOS_FORCE_INLINE
real64
TableFunction::KernelWrapper::getCoord( IN_ARRAY const & input, localIndex const dim, InterpolationType interpolationMethod ) const
{
  // Determine the index to the nearest table entry
  localIndex subIndex;
  arraySlice1d< real64 const > const coords = m_coordinates[dim];
  // Determine the index along each table axis
  if( input[dim] <= coords[0] )
  {
    // Coordinate is to the left of the table axis
    subIndex = 0;
  }
  else if( input[dim] >= coords[coords.size() - 1] )
  {
    // Coordinate is to the right of the table axis
    subIndex = coords.size() - 1;
  }
  else
  {
    // Coordinate is within the table axis
    // Note: find() will return the index of the upper table vertex
    auto const lower = LvArray::sortedArrayManipulation::find( coords.begin(), coords.size(), input[dim] );
    subIndex = LvArray::integerConversion< localIndex >( lower );
 
    // Interpolation types:
    //   - Nearest returns the value of the closest table vertex
    //   - Upper returns the value of the next table vertex
    //   - Lower returns the value of the previous table vertex
    if( interpolationMethod == TableFunction::InterpolationType::Nearest )
    {
      if( ( input[dim] - coords[subIndex - 1]) <= ( coords[subIndex] - input[dim]) )
      {
        --subIndex;
      }
    }
    else if( interpolationMethod == TableFunction::InterpolationType::Lower )
    {
      if( subIndex > 0 )
      {
        --subIndex;
      }
    }
  }
 
  // Retrieve the nearest coordinate
  return coords[subIndex];
}
 
template< typename IN_ARRAY, typename OUT_ARRAY >
GEOS_HOST_DEVICE
GEOS_FORCE_INLINE
real64
TableFunction::KernelWrapper::compute( IN_ARRAY const & input, OUT_ARRAY && derivatives ) const
{
  // Linear interpolation
  if( m_interpolationMethod == TableFunction::InterpolationType::Linear )
  {
    return interpolateLinear( input, derivatives );
  }
  // Nearest, Upper, Lower interpolation methods
  else
  {
    return interpolateRound( input, derivatives );
  }
}
 
template< typename IN_ARRAY, typename OUT_ARRAY >
GEOS_HOST_DEVICE
GEOS_FORCE_INLINE
real64
TableFunction::KernelWrapper::interpolateLinear( IN_ARRAY const & input, OUT_ARRAY && derivatives ) const
{
  integer const numDimensions = LvArray::integerConversion< integer >( m_coordinates.size() );
 
  localIndex bounds[maxDimensions][2]{};
  real64 weights[maxDimensions][2]{};
  real64 dWeights_dInput[maxDimensions][2]{};
 
  // Determine position, weights
  for( integer dim = 0; dim < numDimensions; ++dim )
  {
    arraySlice1d< real64 const > const coords = m_coordinates[dim];
    if( input[dim] <= coords[0] )
    {
      // Coordinate is to the left of this axis
      bounds[dim][0] = 0;
      bounds[dim][1] = 0;
      weights[dim][0] = 0;
      weights[dim][1] = 1;
      dWeights_dInput[dim][0] = 0;
      dWeights_dInput[dim][1] = 0;
    }
    else if( input[dim] >= coords[coords.size() - 1] )
    {
      // Coordinate is to the right of this axis
      bounds[dim][0] = coords.size() - 1;
      bounds[dim][1] = bounds[dim][0];
      weights[dim][0] = 1;
      weights[dim][1] = 0;
      dWeights_dInput[dim][0] = 0;
      dWeights_dInput[dim][1] = 0;
    }
    else
    {
      // Find the coordinate index
      // Note: find uses a binary search...  If we assume coordinates are
      // evenly spaced, we can speed things up considerably
      auto lower = LvArray::sortedArrayManipulation::find( coords.begin(), coords.size(), input[dim] );
      bounds[dim][1] = LvArray::integerConversion< localIndex >( lower );
      bounds[dim][0] = bounds[dim][1] - 1;
 
      real64 const dx = coords[bounds[dim][1]] - coords[bounds[dim][0]];
      weights[dim][0] = 1.0 - ( input[dim] - coords[bounds[dim][0]]) / dx;
      weights[dim][1] = 1.0 - weights[dim][0];
      dWeights_dInput[dim][0] = -1.0 / dx;
      dWeights_dInput[dim][1] = -dWeights_dInput[dim][0];
    }
  }
 
  // Calculate the result
  real64 value = 0.0;
  for( integer dim = 0; dim < numDimensions; ++dim )
  {
    derivatives[dim] = 0.0;
  }
 
  integer const numCorners = 1 << numDimensions;
  for( integer point = 0; point < numCorners; ++point )
  {
    // Find array index
    localIndex tableIndex = 0;
    localIndex stride = 1;
    for( integer dim = 0; dim < numDimensions; ++dim )
    {
      integer const corner = (point >> dim) & 1;
      tableIndex += bounds[dim][corner] * stride;
      stride *= m_coordinates.sizeOfArray( dim );
    }
 
    // Determine weighted value
    real64 cornerValue = m_values[tableIndex];
    real64 dCornerValue_dInput[maxDimensions]{};
    for( integer dim = 0; dim < numDimensions; ++dim )
    {
      dCornerValue_dInput[dim] = cornerValue;
    }
 
    for( integer dim = 0; dim < numDimensions; ++dim )
    {
      integer const corner = (point >> dim) & 1;
      cornerValue *= weights[dim][corner];
      for( integer kk = 0; kk < numDimensions; ++kk )
      {
        dCornerValue_dInput[kk] *= ( dim == kk ) ? dWeights_dInput[dim][corner] : weights[dim][corner];
      }
    }
 
    for( integer dim = 0; dim < numDimensions; ++dim )
    {
      derivatives[dim] += dCornerValue_dInput[dim];
    }
    value += cornerValue;
  }
  return value;
}
 
template< typename IN_ARRAY, typename OUT_ARRAY >
GEOS_HOST_DEVICE
GEOS_FORCE_INLINE
real64
TableFunction::KernelWrapper::interpolateRound( IN_ARRAY const & input, OUT_ARRAY && derivatives ) const
{
  GEOS_UNUSED_VAR( input, derivatives );
  GEOS_ERROR( "Rounding interpolation with derivatives not implemented" );
  return 0.0;
}
 
 
ENUM_STRINGS( TableFunction::InterpolationType,
              "linear",
              "nearest",
              "upper",
              "lower" );
 
template<>
string TableTextFormatter::toString< TableFunction >( TableFunction const & tableData ) const;
 
template<>
string TableCSVFormatter::toString< TableFunction >( TableFunction const & tableData ) const;
 
} /* namespace geos */
 
#endif /* GEOS_FUNCTIONS_TABLEFUNCTION_HPP_ */