FunctionBase.hpp

Go to the documentation of this file.

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2018-2020 Total, S.A
 * Copyright (c) 2019-     GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */

#ifndef GEOSX_MANAGERS_FUNCTIONS_FUNCTIONBASE_HPP_
#define GEOSX_MANAGERS_FUNCTIONS_FUNCTIONBASE_HPP_

#include "dataRepository/Group.hpp"
#include "rajaInterface/GEOS_RAJA_Interface.hpp"

namespace geosx
{

namespace dataRepository
{
namespace keys
{
string const inputVarNames( "inputVarNames" );
}
}

class FunctionBase : public dataRepository::Group
{
public:
  FunctionBase( const std::string & name,
                dataRepository::Group * const parent );

  virtual ~FunctionBase() override;

  static string CatalogName() { return "FunctionBase"; }

  virtual void InitializeFunction(){}

  integer isFunctionOfTime() const;

  virtual void Evaluate( dataRepository::Group const * const group,
                         real64 const time,
                         SortedArrayView< localIndex const > const & set,
                         real64_array & result ) const = 0;

  virtual real64 Evaluate( real64 const * const input ) const = 0;

  using CatalogInterface = dataRepository::CatalogInterface< FunctionBase, std::string const &, Group * const >;

  static CatalogInterface::CatalogType & GetCatalog()
  {
    static CatalogInterface::CatalogType catalog;
    return catalog;
  }

  real64_array EvaluateStats( dataRepository::Group const * const group,
                              real64 const time,
                              SortedArray< localIndex > const & set ) const;

  void setInputVarNames( string_array inputVarNames ) { m_inputVarNames = inputVarNames; }


protected:
  string_array m_inputVarNames;

  template< typename LEAF >
  void EvaluateT( dataRepository::Group const * const group,
                  real64 const time,
                  SortedArrayView< localIndex const > const & set,
                  real64_array & result ) const;

  virtual void PostProcessInput() override { InitializeFunction(); }

};

template< typename LEAF >
void FunctionBase::EvaluateT( dataRepository::Group const * const group,
                              real64 const time,
                              SortedArrayView< localIndex const > const & set,
                              real64_array & result ) const
{
  real64 const * input_ptrs[4];
  localIndex varSize[4] = {0, 0, 0, 0};
  int timeVar[4] = {1, 1, 1, 1};

  arrayView1d< string const > const & inputVarNames = this->getReference< string_array >( dataRepository::keys::inputVarNames );
  localIndex const numVars = LvArray::integerConversion< localIndex >( inputVarNames.size());
  localIndex groupSize = group->size();
  localIndex totalVarSize = 0;
  for( auto varIndex=0; varIndex<numVars; ++varIndex )
  {
    string const & varName = inputVarNames[varIndex];

    if( varName=="time" )
    {
      input_ptrs[varIndex] = &time;
      varSize[varIndex] = 1;
      timeVar[varIndex] = 0;
      ++totalVarSize;
    }
    else if( groupSize > 0 )
    {
      // Should we throw a warning if the group is zero-length?
      dataRepository::WrapperBase const * wrapper = group->getWrapperBase( varName );
      input_ptrs[ varIndex ] = reinterpret_cast< double const * >( wrapper->voidPointer() );

      localIndex wrapperSize = LvArray::integerConversion< localIndex >( wrapper->size());
      varSize[varIndex] = wrapperSize / groupSize;
      totalVarSize += varSize[varIndex];
    }
  }

  // Make sure the inputs do not exceed the maximum length
  GEOSX_ERROR_IF( totalVarSize > 4, "Function input size is: " << totalVarSize );

  // Make sure the result / set size match
  GEOSX_ERROR_IF( result.size() != set.size(), "To apply a function to a set, the size of the result and set must match" );


  forAll< serialPolicy >( set.size(), [&, set]( localIndex const i )
  {
    localIndex const index = set[ i ];
    double input[4];
    int c = 0;
    for( int a=0; a<numVars; ++a )
    {
      for( int b=0; b<varSize[a]; ++b )
      {
        input[c] = input_ptrs[a][(index*varSize[a]+b)*timeVar[a]];
        ++c;
      }
    }

    // Note: we expect that result is the same size as the set
    result[i] = static_cast< LEAF const * >(this)->Evaluate( input );
  } );
}
} /* namespace geosx */

#endif /* GEOSX_MANAGERS_FUNCTIONS_FUNCTIONBASE_HPP_ */