FieldSpecificationBase.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_FIELDSPECIFICATION_FIELDSPECIFICATIONBASE_HPP
#define GEOS_FIELDSPECIFICATION_FIELDSPECIFICATIONBASE_HPP
 
#include "common/DataTypes.hpp"
#include "common/TypeDispatch.hpp"
#include "codingUtilities/traits.hpp"
#include "codingUtilities/Utilities.hpp"
#include "dataRepository/Group.hpp"
#include "functions/FunctionBase.hpp"
#include "common/FieldSpecificationOps.hpp"
#include "mesh/ObjectManagerBase.hpp"
#include "mesh/MeshObjectPath.hpp"
#include "functions/FunctionManager.hpp"
#include "common/GEOS_RAJA_Interface.hpp"
 
namespace geos
{
class Function;
 
 
class FieldSpecificationBase : public dataRepository::Group
{
public:
 
  using CatalogInterface = dataRepository::CatalogInterface< FieldSpecificationBase,
                                                             string const &,
                                                             dataRepository::Group * const >;
  enum class SetErrorMode : integer
  {
    silent,
    error,
    warning
  };
 
  static CatalogInterface::CatalogType & getCatalog();
 
  static string catalogName() { return "FieldSpecification"; }
 
  virtual const string getCatalogName() const
  {
    return FieldSpecificationBase::catalogName();
  }
 
  FieldSpecificationBase( string const & name, dataRepository::Group * parent );
 
  virtual ~FieldSpecificationBase() override;
 
 
  FieldSpecificationBase( FieldSpecificationBase const & ) = delete;
 
  FieldSpecificationBase( FieldSpecificationBase && ) = default;
 
  FieldSpecificationBase & operator=( FieldSpecificationBase const & ) = delete;
 
  FieldSpecificationBase & operator=( FieldSpecificationBase && ) = delete;
 
  template< typename OBJECT_TYPE,
            typename BC_TYPE = FieldSpecificationBase,
            typename LAMBDA >
  void apply( MeshLevel & mesh,
              LAMBDA && lambda ) const
  {
    MeshObjectPath const & meshObjectPaths = this->getMeshObjectPaths();
    meshObjectPaths.forObjectsInPath< OBJECT_TYPE >( mesh,
                                                     [&] ( OBJECT_TYPE & object )
    {
      {
        dataRepository::Group const & setGroup = object.getGroup( ObjectManagerBase::groupKeyStruct::setsString() );
        string_array setNames = this->getSetNames();
        for( auto & setName : setNames )
        {
          if( setGroup.hasWrapper( setName ) )
          {
            SortedArrayView< localIndex const > const & targetSet = setGroup.getReference< SortedArray< localIndex > >( setName );
            lambda( dynamic_cast< BC_TYPE const & >(*this), setName, targetSet, object, getFieldName() );
          }
        }
      }
    } );
  }
 
  template< typename FIELD_OP, typename POLICY, typename T, int N, int USD >
  void applyFieldValueKernel( ArrayView< T, N, USD > const & field,
                              SortedArrayView< localIndex const > const & targetSet,
                              real64 const time,
                              Group & dataGroup ) const;
 
  template< typename FIELD_OP, typename POLICY=parallelHostPolicy >
  void applyFieldValue( SortedArrayView< localIndex const > const & targetSet,
                        real64 const time,
                        dataRepository::Group & dataGroup,
                        string const & fieldname ) const;
 
  template< typename FIELD_OP, typename POLICY, typename T, int NDIM, int USD >
  void applyBoundaryConditionToSystemKernel( SortedArrayView< localIndex const > const & targetSet,
                                             real64 const time,
                                             dataRepository::Group const & dataGroup,
                                             arrayView1d< globalIndex const > const & dofMap,
                                             globalIndex const dofRankOffset,
                                             CRSMatrixView< real64, globalIndex const > const & matrix,
                                             arrayView1d< real64 > const & rhs,
                                             ArrayView< T const, NDIM, USD > const & fieldView ) const;
 
  template< typename FIELD_OP, typename POLICY >
  void applyBoundaryConditionToSystem( SortedArrayView< localIndex const > const & targetSet,
                                       real64 const time,
                                       dataRepository::Group const & dataGroup,
                                       string const & fieldName,
                                       string const & dofMapName,
                                       globalIndex const dofRankOffset,
                                       CRSMatrixView< real64, globalIndex const > const & matrix,
                                       arrayView1d< real64 > const & rhs ) const;
 
  template< typename FIELD_OP, typename POLICY, typename LAMBDA >
  void
  applyBoundaryConditionToSystem( SortedArrayView< localIndex const > const & targetSet,
                                  real64 const time,
                                  dataRepository::Group const & dataGroup,
                                  arrayView1d< globalIndex const > const & dofMap,
                                  globalIndex const dofRankOffset,
                                  CRSMatrixView< real64, globalIndex const > const & matrix,
                                  arrayView1d< real64 > const & rhs,
                                  LAMBDA && lambda ) const;
 
  template< typename FIELD_OP, typename POLICY, typename LAMBDA >
  void
  applyBoundaryConditionToSystem( SortedArrayView< localIndex const > const & targetSet,
                                  real64 const time,
                                  real64 const dt,
                                  dataRepository::Group const & dataGroup,
                                  arrayView1d< globalIndex const > const & dofMap,
                                  globalIndex const dofRankOffset,
                                  CRSMatrixView< real64, globalIndex const > const & matrix,
                                  arrayView1d< real64 > const & rhs,
                                  LAMBDA && lambda ) const;
 
  template< typename FIELD_OP, typename POLICY, typename LAMBDA >
  void
  computeRhsContribution( SortedArrayView< localIndex const > const & targetSet,
                          real64 const time,
                          real64 const dt,
                          dataRepository::Group const & dataGroup,
                          arrayView1d< globalIndex const > const & dofMap,
                          globalIndex const dofRankOffset,
                          CRSMatrixView< real64, globalIndex const > const & matrix,
                          arrayView1d< globalIndex > const & dof,
                          arrayView1d< real64 > const & rhsContribution,
                          LAMBDA && lambda ) const;
 
 
  template< typename POLICY >
  void zeroSystemRowsForBoundaryCondition( SortedArrayView< localIndex const > const & targetSet,
                                           arrayView1d< globalIndex const > const & dofMap,
                                           CRSMatrixView< real64, globalIndex const > const & matrix ) const;
 
  struct viewKeyStruct
  {
    constexpr static char const * setNamesString() { return "setNames"; }
    constexpr static char const * constitutivePathString() { return "constitutivePath"; }
    constexpr static char const * objectPathString() { return "objectPath"; }
    constexpr static char const * fieldNameString() { return "fieldName"; }
    constexpr static char const * dataTypeString() { return "dataType"; }
    constexpr static char const * componentString() { return "component"; }
    constexpr static char const * directionString() { return "direction"; }
    constexpr static char const * bcApplicationTableNameString() { return "bcApplicationTableName"; }
    constexpr static char const * scaleString() { return "scale"; }
    constexpr static char const * functionNameString() { return "functionName"; }
    constexpr static char const * initialConditionString() { return "initialCondition"; }
    constexpr static char const * beginTimeString() { return "beginTime"; }
    constexpr static char const * endTimeString() { return "endTime"; }
    constexpr static char const * errorSetModeString() { return "errorSetMode"; }
  };
 
  string const & getFunctionName() const
  {
    return m_functionName;
  }
 
  virtual const string & getObjectPath() const
  {
    return m_objectPath;
  }
 
  virtual const string & getFieldName() const
  {
    return m_fieldName;
  }
 
  virtual int getComponent() const
  {
    return m_component;
  }
 
  virtual R1Tensor const & getDirection() const
  {
    GEOS_UNUSED_VAR( time );
    return m_direction;
  }
 
  real64 getStartTime() const
  {
    return m_beginTime;
  }
 
  real64 getEndTime() const
  {
    return m_endTime;
  }
 
  string_array const & getSetNames() const
  {
    return m_setNames;
  }
 
  int initialCondition() const
  {
    return m_initialCondition;
  }
 
  real64 getScale() const
  {
    return m_scale;
  }
 
  void setFieldName( string const & fieldName )
  {
    m_fieldName = fieldName;
  }
 
  void setObjectPath( string const & objectPath )
  {
    m_objectPath = objectPath;
  }
 
  void setScale( real64 const & scale )
  {
    m_scale = scale;
  }
 
  void initialCondition( bool isInitialCondition )
  {
    m_initialCondition = isInitialCondition;
  }
 
  void addSetName( string const & setName )
  {
    m_setNames.emplace_back( setName );
  }
 
  void setMeshObjectPath( Group const & meshBodies );
 
  MeshObjectPath const & getMeshObjectPaths() const
  {
    return *(m_meshObjectPaths.get());
  }
 
 
protected:
 
 
private:
 
 
  string_array m_setNames;
 
  string m_objectPath;
 
  std::unique_ptr< MeshObjectPath > m_meshObjectPaths;
 
  string m_fieldName;
 
 
  int m_component;
 
  R1Tensor m_direction;
 
  int m_initialCondition;
 
  string m_functionName;
 
  real64 m_scale;
 
  real64 m_beginTime;
 
  real64 m_endTime;
 
  string m_bcApplicationFunctionName;
 
  SetErrorMode m_emptySetErrorMode;
};
 
 
template< typename FIELD_OP, typename POLICY, typename T, int N, int USD >
void FieldSpecificationBase::applyFieldValueKernel( ArrayView< T, N, USD > const & field,
                                                    SortedArrayView< localIndex const > const & targetSet,
                                                    real64 const time,
                                                    Group & dataGroup ) const
{
  integer const component = getComponent();
  FunctionManager & functionManager = FunctionManager::getInstance();
 
  if( m_functionName.empty() )
  {
    real64 const value = m_scale;
    forAll< POLICY >( targetSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const i )
    {
      localIndex const a = targetSet[ i ];
      FIELD_OP::SpecifyFieldValue( field, a, component, value );
    } );
  }
  else
  {
    FunctionBase const & function = [&]() -> FunctionBase const &
    {
      try
      {
        return functionManager.getGroup< FunctionBase >( m_functionName );
      }
      catch( std::exception const & e )
      {
        throw InputError( e, GEOS_FMT( "Error while reading {}:\n",
                                       getWrapperDataContext( viewKeyStruct::functionNameString() ) ) );
      }
    }();
 
    if( function.isFunctionOfTime()==2 )
    {
      real64 const value = m_scale * function.evaluate( &time );
      forAll< POLICY >( targetSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const i )
      {
        localIndex const a = targetSet[ i ];
        FIELD_OP::SpecifyFieldValue( field, a, component, value );
      } );
    }
    else
    {
      real64_array result( static_cast< localIndex >( targetSet.size() ) );
      function.evaluate( dataGroup, time, targetSet, result );
      arrayView1d< real64 const > const & resultView = result.toViewConst();
      real64 const scale = m_scale;
      forAll< POLICY >( targetSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const i )
      {
        localIndex const a = targetSet[ i ];
        FIELD_OP::SpecifyFieldValue( field, a, component, scale * resultView[i] );
      } );
    }
  }
}
 
 
template< typename FIELD_OP, typename POLICY >
void FieldSpecificationBase::applyFieldValue( SortedArrayView< localIndex const > const & targetSet,
                                              real64 const time,
                                              dataRepository::Group & dataGroup,
                                              string const & fieldName ) const
{
  dataRepository::WrapperBase & wrapper = dataGroup.getWrapperBase( fieldName );
 
  // // This function is used in setting boundary/initial conditions on simulation fields.
  // // This is meaningful for 1/2/3D real arrays and sometimes 1D integer (indicator) arrays.
  using FieldTypes = types::ListofTypeList< types::Join< types::ArrayTypes< types::RealTypes, types::DimsUpTo< 3 > >,
                                                         types::ArrayTypes< types::TypeList< integer >, types::DimsSingle< 1 > > > >;
 
 
  types::dispatch( FieldTypes{}, [&]( auto tupleOfTypes )
  {
    using ArrayType = camp::first< decltype( tupleOfTypes ) >;
    auto & wrapperT = dataRepository::Wrapper< ArrayType >::cast( wrapper );
    applyFieldValueKernel< FIELD_OP, POLICY >( wrapperT.reference().toView(), targetSet, time, dataGroup );
  }, wrapper );
}
 
template< typename FIELD_OP, typename POLICY, typename T, int NDIM, int USD >
void FieldSpecificationBase::applyBoundaryConditionToSystemKernel( SortedArrayView< localIndex const > const & targetSet,
                                                                   real64 const time,
                                                                   dataRepository::Group const & dataGroup,
                                                                   arrayView1d< globalIndex const > const & dofMap,
                                                                   globalIndex const dofRankOffset,
                                                                   CRSMatrixView< real64, globalIndex const > const & matrix,
                                                                   arrayView1d< real64 > const & rhs,
                                                                   ArrayView< T const, NDIM, USD > const & fieldView ) const
{
  integer const component = getComponent();
  this->applyBoundaryConditionToSystem< FIELD_OP, POLICY >( targetSet, time, dataGroup, dofMap, dofRankOffset, matrix, rhs,
                                                            [fieldView, component] GEOS_HOST_DEVICE ( localIndex const a )
  {
    real64 value = 0.0;
    FieldSpecificationEqual::readFieldValue( fieldView, a, component, value );
    return value;
  } );
}
 
template< typename FIELD_OP, typename POLICY >
void FieldSpecificationBase::applyBoundaryConditionToSystem( SortedArrayView< localIndex const > const & targetSet,
                                                             real64 const time,
                                                             dataRepository::Group const & dataGroup,
                                                             string const & fieldName,
                                                             string const & dofMapName,
                                                             globalIndex const dofRankOffset,
                                                             CRSMatrixView< real64, globalIndex const > const & matrix,
                                                             arrayView1d< real64 > const & rhs ) const
{
  dataRepository::WrapperBase const & wrapper = dataGroup.getWrapperBase( fieldName );
  arrayView1d< globalIndex const > const & dofMap = dataGroup.getReference< array1d< globalIndex > >( dofMapName );
 
  // We're reading values from a field, which is only well-defined for dims 1 and 2
  using FieldTypes = types::ListofTypeList< types::ArrayTypes< types::RealTypes, types::DimsUpTo< 2 > > >;
  types::dispatch( FieldTypes{}, [&]( auto tupleOfTypes )
  {
    using ArrayType = camp::first< decltype( tupleOfTypes ) >;
    auto const & wrapperT = dataRepository::Wrapper< ArrayType >::cast( wrapper );
    applyBoundaryConditionToSystemKernel< FIELD_OP, POLICY >( targetSet,
                                                              time,
                                                              dataGroup,
                                                              dofMap,
                                                              dofRankOffset,
                                                              matrix,
                                                              rhs,
                                                              wrapperT.reference() );
  }, wrapper );
}
 
template< typename FIELD_OP, typename POLICY, typename LAMBDA >
void
FieldSpecificationBase::
  applyBoundaryConditionToSystem( SortedArrayView< localIndex const > const & targetSet,
                                  real64 const time,
                                  dataRepository::Group const & dataGroup,
                                  arrayView1d< globalIndex const > const & dofMap,
                                  globalIndex const dofRankOffset,
                                  CRSMatrixView< real64, globalIndex const > const & matrix,
                                  arrayView1d< real64 > const & rhs,
                                  LAMBDA && lambda ) const
{
  return applyBoundaryConditionToSystem< FIELD_OP, POLICY >( targetSet,
                                                             time,
                                                             1.0,
                                                             dataGroup,
                                                             dofMap,
                                                             dofRankOffset,
                                                             matrix,
                                                             rhs,
                                                             std::forward< LAMBDA >( lambda ) );
}
 
template< typename FIELD_OP, typename POLICY, typename LAMBDA >
void
FieldSpecificationBase::
  applyBoundaryConditionToSystem( SortedArrayView< localIndex const > const & targetSet,
                                  real64 const time,
                                  real64 const dt,
                                  dataRepository::Group const & dataGroup,
                                  arrayView1d< globalIndex const > const & dofMap,
                                  globalIndex const dofRankOffset,
                                  CRSMatrixView< real64, globalIndex const > const & matrix,
                                  arrayView1d< real64 > const & rhs,
                                  LAMBDA && lambda ) const
{
  array1d< globalIndex > dofArray( targetSet.size() );
  arrayView1d< globalIndex > const & dof = dofArray.toView();
 
  array1d< real64 > rhsContributionArray( targetSet.size() );
  arrayView1d< real64 > const & rhsContribution = rhsContributionArray.toView();
 
  computeRhsContribution< FIELD_OP, POLICY, LAMBDA >( targetSet,
                                                      time,
                                                      dt,
                                                      dataGroup,
                                                      dofMap,
                                                      dofRankOffset,
                                                      matrix,
                                                      dof,
                                                      rhsContribution,
                                                      std::forward< LAMBDA >( lambda ) );
 
  FIELD_OP::template prescribeRhsValues< POLICY >( rhs, dof, dofRankOffset, rhsContribution );
}
 
template< typename FIELD_OP, typename POLICY, typename LAMBDA >
void
FieldSpecificationBase::
  computeRhsContribution( SortedArrayView< localIndex const > const & targetSet,
                          real64 const time,
                          real64 const dt,
                          dataRepository::Group const & dataGroup,
                          arrayView1d< globalIndex const > const & dofMap,
                          globalIndex const dofRankOffset,
                          CRSMatrixView< real64, globalIndex const > const & matrix,
                          arrayView1d< globalIndex > const & dof,
                          arrayView1d< real64 > const & rhsContribution,
                          LAMBDA && lambda ) const
{
  integer const component = ( getComponent() >=0 ) ? getComponent() : 0;
  string const & functionName = getReference< string >( viewKeyStruct::functionNameString() );
  FunctionManager & functionManager = FunctionManager::getInstance();
 
  // Compute the value of the rhs terms, and collect the dof numbers
  // The rhs terms will be assembled in applyBoundaryConditionToSystem (or in the solver for CompositionalMultiphaseBase)
 
  if( functionName.empty() || functionManager.getGroup< FunctionBase >( functionName ).isFunctionOfTime() == 2 )
  {
    real64 value = m_scale * dt;
    if( !functionName.empty() )
    {
      FunctionBase const & function = functionManager.getGroup< FunctionBase >( functionName );
      value *= function.evaluate( &time );
    }
 
    forAll< POLICY >( targetSet.size(),
                      [targetSet, dof, dofMap, dofRankOffset, component, matrix, rhsContribution, value, lambda] GEOS_HOST_DEVICE ( localIndex const i )
    {
      localIndex const a = targetSet[ i ];
      dof[ i ] = dofMap[ a ] + component;
      FIELD_OP::SpecifyFieldValue( dof[ i ],
                                   dofRankOffset,
                                   matrix,
                                   rhsContribution[ i ],
                                   value,
                                   lambda( a ) );
    } );
  }
  else
  {
    FunctionBase const & function = functionManager.getGroup< FunctionBase >( functionName );
 
    real64_array resultsArray( targetSet.size() );
    function.evaluate( dataGroup, time, targetSet, resultsArray );
    arrayView1d< real64 const > const & results = resultsArray.toViewConst();
    real64 const value = m_scale * dt;
 
    forAll< POLICY >( targetSet.size(),
                      [targetSet, dof, dofMap, dofRankOffset, component, matrix, rhsContribution, results, value, lambda] GEOS_HOST_DEVICE (
                        localIndex const i )
    {
      localIndex const a = targetSet[ i ];
      dof[ i ] = dofMap[ a ] + component;
      FIELD_OP::SpecifyFieldValue( dof[ i ],
                                   dofRankOffset,
                                   matrix,
                                   rhsContribution[ i ],
                                   value * results[ i ],
                                   lambda( a ) );
    } );
  }
}
 
 
template< typename POLICY >
void FieldSpecificationBase::zeroSystemRowsForBoundaryCondition( SortedArrayView< localIndex const > const & targetSet,
                                                                 arrayView1d< globalIndex const > const & dofMap,
                                                                 CRSMatrixView< real64, globalIndex const > const & matrix ) const
 
{
  integer const component = ( getComponent() >=0 ) ? getComponent() : 0;
  forAll< POLICY >( targetSet.size(), [targetSet, dofMap, matrix, component] GEOS_HOST_DEVICE ( localIndex const i )
  {
    localIndex const a = targetSet[ i ];
    globalIndex const dof = dofMap[ a ] + component;
 
    arraySlice1d< real64 > const entries = matrix.getEntries( dof );
    localIndex const numEntries = matrix.numNonZeros( dof );
 
    for( localIndex j = 0; j < numEntries; ++j )
    {
      entries[ j ] = 0;
    }
  } );
}
 
 
ENUM_STRINGS( FieldSpecificationBase::SetErrorMode,
              "silent",
              "error",
              "warning" );
 
}
 
#endif //GEOS_FIELDSPECIFICATION_FIELDSPECIFICATIONBASE_HPP