FieldSpecificationBase.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_FIELDSPECIFICATION_FIELDSPECIFICATIONBASE_HPP
#define GEOSX_MANAGERS_FIELDSPECIFICATION_FIELDSPECIFICATIONBASE_HPP

#include "common/DataTypes.hpp"
#include "codingUtilities/traits.hpp"
#include "codingUtilities/Utilities.hpp"
#include "dataRepository/Group.hpp"
#include "linearAlgebra/interfaces/InterfaceTypes.hpp"
#include "managers/FieldSpecification/FieldSpecificationOps.hpp"
#include "managers/Functions/FunctionManager.hpp"
#include "rajaInterface/GEOS_RAJA_Interface.hpp"
#include "managers/ObjectManagerBase.hpp"


namespace geosx
{
class Function;


class FieldSpecificationBase : public dataRepository::Group
{
public:

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

  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;

  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 LAI >
  void ApplyBoundaryConditionToSystem( SortedArrayView< localIndex const > const & targetSet,
                                       real64 const time,
                                       dataRepository::Group const * const dataGroup,
                                       string const & fieldName,
                                       string const & dofMapName,
                                       integer const & dofDim,
                                       typename LAI::ParallelMatrix & matrix,
                                       typename LAI::ParallelVector & rhs ) const;


  template< typename FIELD_OP, typename LAI, typename LAMBDA >
  void
  ApplyBoundaryConditionToSystem( SortedArrayView< localIndex const > const & targetSet,
                                  real64 const time,
                                  dataRepository::Group const * const dataGroup,
                                  arrayView1d< globalIndex const > const & dofMap,
                                  integer const & dofDim,
                                  typename LAI::ParallelMatrix & matrix,
                                  typename LAI::ParallelVector & rhs,
                                  LAMBDA && lambda ) const;

  template< typename FIELD_OP, typename LAI, typename LAMBDA >
  void
  ApplyBoundaryConditionToSystem( SortedArrayView< localIndex const > const & targetSet,
                                  real64 const time,
                                  real64 const dt,
                                  dataRepository::Group const * const dataGroup,
                                  arrayView1d< globalIndex const > const & dofMap,
                                  integer const & dofDim,
                                  typename LAI::ParallelMatrix & matrix,
                                  typename LAI::ParallelVector & rhs,
                                  LAMBDA && lambda ) const;

  template< typename LAI >
  void ZeroSystemRowsForBoundaryCondition( SortedArrayView< localIndex const > const & targetSet,
                                           arrayView1d< globalIndex const > const & dofMap,
                                           typename LAI::ParallelMatrix & matrix ) 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 * 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 * 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 * 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 * 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 POLICY >
  void ZeroSystemRowsForBoundaryCondition( SortedArrayView< localIndex const > const & targetSet,
                                           arrayView1d< globalIndex const > const & dofMap,
                                           CRSMatrixView< real64, globalIndex const > const & matrix ) const;

  struct viewKeyStruct
  {
    constexpr static auto setNamesString = "setNames";
    constexpr static auto constitutivePathString = "constitutivePath";
    constexpr static auto objectPathString = "objectPath";
    constexpr static auto fieldNameString = "fieldName";
    constexpr static auto dataTypeString = "dataType";
    constexpr static auto componentString = "component";
    constexpr static auto directionString = "direction";
    constexpr static auto bcApplicationTableNameString = "bcApplicationTableName";
    constexpr static auto scaleString = "scale";
    constexpr static auto functionNameString = "functionName";
    constexpr static auto initialConditionString = "initialCondition";
    constexpr static auto beginTimeString = "beginTime";
    constexpr static auto endTimeString = "endTime";
    constexpr static auto fluxBoundaryConditionString = "fluxBoundaryCondition";
  } viewKeys; 

  struct groupKeyStruct
  {} groupKeys; 


  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 const R1Tensor & GetDirection( real64 time )
  {
    GEOSX_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 );
  }


protected:
  void PostProcessInput() override final;

  bool m_normalizeBySetSize;

private:


  string_array m_setNames;

  string m_objectPath;

  string m_fieldName;


//  string m_dataType;

  int m_component;

  R1Tensor m_direction;

  int m_initialCondition;

  string m_functionName;

  real64 m_scale;

  real64 m_beginTime;

  real64 m_endTime;

  string m_bcApplicationFunctionName;

  //real64 m_setSizeScalingFactor;

};


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::Instance();

  if( m_functionName.empty() )
  {
    real64 const value = m_scale;
    forAll< POLICY >( targetSet.size(), [=] GEOSX_HOST_DEVICE ( localIndex const i )
    {
      localIndex const a = targetSet[ i ];
      FIELD_OP::SpecifyFieldValue( field, a, component, value );
    } );
  }
  else
  {
    FunctionBase const * const function  = functionManager.GetGroup< FunctionBase >( m_functionName );

    GEOSX_ERROR_IF( function == nullptr, "Function '" << m_functionName << "' not found" );

    if( function->isFunctionOfTime()==2 )
    {
      real64 const value = m_scale * function->Evaluate( &time );
      forAll< POLICY >( targetSet.size(), [=] GEOSX_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(), [=] GEOSX_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 );
  std::type_index typeIndex = std::type_index( wrapper->get_typeid());

  rtTypes::ApplyArrayTypeLambda2( rtTypes::typeID( typeIndex ),
                                  true,
                                  [&]( auto arrayInstance, auto GEOSX_UNUSED_PARAM( dataTypeInstance ) )
  {
    using ArrayType = decltype(arrayInstance);
    dataRepository::Wrapper< ArrayType > & view = dataRepository::Wrapper< ArrayType >::cast( *wrapper );

    auto const & field = view.reference().toView();
    ApplyFieldValueKernel< FIELD_OP, POLICY >( field, targetSet, time, dataGroup );
  } );
}

template< typename FIELD_OP, typename LAI >
void FieldSpecificationBase::ApplyBoundaryConditionToSystem( SortedArrayView< localIndex const > const & targetSet,
                                                             real64 const time,
                                                             dataRepository::Group const * const dataGroup,
                                                             string const & fieldName,
                                                             string const & dofMapName,
                                                             integer const & dofDim,
                                                             typename LAI::ParallelMatrix & matrix,
                                                             typename LAI::ParallelVector & rhs ) const
{
  dataRepository::WrapperBase const & wrapperBase = *dataGroup->getWrapperBase( fieldName );
  std::type_index typeIndex = std::type_index( wrapperBase.get_typeid());
  arrayView1d< globalIndex const > const & dofMap = dataGroup->getReference< array1d< globalIndex > >( dofMapName );
  integer const component = GetComponent();

  rtTypes::ApplyArrayTypeLambda1( rtTypes::typeID( typeIndex ), [&]( auto type )
  {
    using FieldType = decltype( type );
    dataRepository::Wrapper< FieldType > const & wrapper = dataRepository::Wrapper< FieldType >::cast( wrapperBase );
    traits::ViewTypeConst< FieldType > fieldView = wrapper.reference();

    this->ApplyBoundaryConditionToSystem< FIELD_OP, LAI >( targetSet, time, dataGroup, dofMap, dofDim, matrix, rhs,
                                                           [&]( localIndex const a )
    {
      real64 value = 0.0;
      FieldSpecificationEqual::ReadFieldValue( fieldView, a, component, value );
      return value;
    } );
  } );
}

template< typename FIELD_OP, typename LAI, typename LAMBDA >
void
FieldSpecificationBase::
  ApplyBoundaryConditionToSystem( SortedArrayView< localIndex const > const & targetSet,
                                  real64 const time,
                                  dataRepository::Group const * const dataGroup,
                                  arrayView1d< globalIndex const > const & dofMap,
                                  integer const & dofDim,
                                  typename LAI::ParallelMatrix & matrix,
                                  typename LAI::ParallelVector & rhs,
                                  LAMBDA && lambda ) const
{
  GEOSX_UNUSED_VAR( dofDim );

  integer const component = GetComponent();
  FunctionManager & functionManager = FunctionManager::Instance();

  globalIndex_array dof( targetSet.size() );
  real64_array rhsContribution( targetSet.size() );

  real64 sizeScalingFactor = 0;
  if( m_normalizeBySetSize )
  {
    // note: this assumes that the ghost elements have been filtered out

    // recompute the set size here to make sure that topology changes are accounted for
    integer const localSetSize = targetSet.size();
    integer globalSetSize = 0;

    // synchronize
    MpiWrapper::allReduce( &localSetSize, &globalSetSize, 1, MPI_SUM, MPI_COMM_GEOSX );

    // set the scaling factor
    sizeScalingFactor = globalSetSize >= 1 ? 1.0 / globalSetSize : 1;
  }
  else
  {
    sizeScalingFactor = 1;
  }

  if( m_functionName.empty() )
  {

    integer counter=0;
    for( auto a : targetSet )
    {
      dof( counter ) = dofMap[a]+component;
      FIELD_OP::template SpecifyFieldValue< LAI >( dof( counter ),
                                                   matrix,
                                                   rhsContribution( counter ),
                                                   m_scale * sizeScalingFactor,
                                                   lambda( a ) );
      ++counter;
    }
    FIELD_OP::template PrescribeRhsValues< LAI >( rhs, counter, dof.data(), rhsContribution.data() );
  }
  else
  {
    FunctionBase const * const function  = functionManager.GetGroup< FunctionBase >( m_functionName );

    GEOSX_ERROR_IF( function == nullptr, "Function '" << m_functionName << "' not found" );

    if( function->isFunctionOfTime()==2 )
    {
      real64 value = m_scale * function->Evaluate( &time ) * sizeScalingFactor;
      integer counter=0;
      for( auto a : targetSet )
      {
        dof( counter ) = dofMap[a] + component;
        FIELD_OP::template SpecifyFieldValue< LAI >( dof( counter ),
                                                     matrix,
                                                     rhsContribution( counter ),
                                                     value,
                                                     lambda( a ) );
        ++counter;
      }
      FIELD_OP::template PrescribeRhsValues< LAI >( rhs, counter, dof.data(), rhsContribution.data() );
    }
    else
    {
      real64_array result;
      result.resize( LvArray::integerConversion< localIndex >( targetSet.size()));
      function->Evaluate( dataGroup, time, targetSet, result );
      integer counter=0;
      for( auto a : targetSet )
      {
        dof( counter ) = dofMap[a] + component;
        FIELD_OP::template SpecifyFieldValue< LAI >( dof( counter ),
                                                     matrix,
                                                     rhsContribution( counter ),
                                                     m_scale * result[counter] * sizeScalingFactor,
                                                     lambda( a ) );
        ++counter;
      }
      FIELD_OP::template PrescribeRhsValues< LAI >( rhs, counter, dof.data(), rhsContribution.data() );
    }
  }
}

template< typename FIELD_OP, typename LAI, typename LAMBDA >
void
FieldSpecificationBase::
  ApplyBoundaryConditionToSystem( SortedArrayView< localIndex const > const & targetSet,
                                  real64 const time,
                                  real64 const dt,
                                  dataRepository::Group const * const dataGroup,
                                  arrayView1d< globalIndex const > const & dofMap,
                                  integer const & dofDim,
                                  typename LAI::ParallelMatrix & matrix,
                                  typename LAI::ParallelVector & rhs,
                                  LAMBDA && lambda ) const
{
  GEOSX_UNUSED_VAR( dofDim );

  integer const component = GetComponent();
  FunctionManager & functionManager = FunctionManager::Instance();

  globalIndex_array dof( targetSet.size() );
  real64_array rhsContribution( targetSet.size() );

  real64 sizeScalingFactor = 0.0;
  if( m_normalizeBySetSize )
  {
    // note: this assumes that the ghost elements have been filtered out

    // recompute the set size here to make sure that topology changes are accounted for
    integer const localSetSize = targetSet.size();
    integer globalSetSize = 0;

    // synchronize
    MpiWrapper::allReduce( &localSetSize, &globalSetSize, 1, MPI_SUM, MPI_COMM_GEOSX );

    // set the scaling factor
    sizeScalingFactor = globalSetSize >= 1 ? 1.0 / globalSetSize : 1;
  }
  else
  {
    sizeScalingFactor = 1;
  }


  if( m_functionName.empty() )
  {

    integer counter=0;
    for( auto a : targetSet )
    {
      dof( counter ) = dofMap[a]+component;
      FIELD_OP::template SpecifyFieldValue< LAI >( dof( counter ),
                                                   matrix,
                                                   rhsContribution( counter ),
                                                   m_scale * dt * sizeScalingFactor,
                                                   lambda( a ) );
      ++counter;
    }
    FIELD_OP::template PrescribeRhsValues< LAI >( rhs, counter, dof.data(), rhsContribution.data() );
  }
  else
  {
    FunctionBase const * const function  = functionManager.GetGroup< FunctionBase >( m_functionName );

    GEOSX_ERROR_IF( function == nullptr, "Function '" << m_functionName << "' not found" );

    if( function->isFunctionOfTime()==2 )
    {
      real64 value = m_scale * dt * function->Evaluate( &time ) * sizeScalingFactor;
      integer counter=0;
      for( auto a : targetSet )
      {
        dof( counter ) = dofMap[a] + component;
        FIELD_OP::template SpecifyFieldValue< LAI >( dof( counter ),
                                                     matrix,
                                                     rhsContribution( counter ),
                                                     value,
                                                     lambda( a ) );
        ++counter;
      }
      FIELD_OP::template PrescribeRhsValues< LAI >( rhs, counter, dof.data(), rhsContribution.data() );
    }
    else
    {
      real64_array result;
      result.resize( LvArray::integerConversion< localIndex >( targetSet.size()));
      function->Evaluate( dataGroup, time, targetSet, result );
      integer counter=0;
      for( auto a : targetSet )
      {
        dof( counter ) = dofMap[a] + component;
        FIELD_OP::template SpecifyFieldValue< LAI >( dof( counter ),
                                                     matrix,
                                                     rhsContribution( counter ),
                                                     m_scale * dt * result[counter] * sizeScalingFactor,
                                                     lambda( a ) );
        ++counter;
      }
      FIELD_OP::template PrescribeRhsValues< LAI >( rhs, counter, dof.data(), rhsContribution.data() );
    }
  }
}

template< typename LAI >
void FieldSpecificationBase::ZeroSystemRowsForBoundaryCondition( SortedArrayView< localIndex const > const & targetSet,
                                                                 arrayView1d< globalIndex const > const & dofMap,
                                                                 typename LAI::ParallelMatrix & matrix ) const

{
  integer const component = GetComponent();
  for( auto a : targetSet )
  {
    globalIndex const dof = dofMap[a]+component;
    matrix.clearRow( dof );
  }
}

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 * 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] GEOSX_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 * 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 & wrapperBase = *dataGroup->getWrapperBase( fieldName );
  std::type_index typeIndex = std::type_index( wrapperBase.get_typeid());
  arrayView1d< globalIndex const > const & dofMap = dataGroup->getReference< array1d< globalIndex > >( dofMapName );

  rtTypes::ApplyArrayTypeLambda1( rtTypes::typeID( typeIndex ), [&]( auto type )
  {
    using FieldType = decltype( type );
    dataRepository::Wrapper< FieldType > const & wrapper = dataRepository::Wrapper< FieldType >::cast( wrapperBase );
    ApplyBoundaryConditionToSystemKernel< FIELD_OP, POLICY >( targetSet, time, dataGroup, dofMap, dofRankOffset, matrix, rhs, wrapper.reference() );
  } );
}

template< typename FIELD_OP, typename POLICY, typename LAMBDA >
void
FieldSpecificationBase::
  ApplyBoundaryConditionToSystem( SortedArrayView< localIndex const > const & targetSet,
                                  real64 const time,
                                  dataRepository::Group const * 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 * const dataGroup,
                                  arrayView1d< globalIndex const > const & dofMap,
                                  globalIndex const dofRankOffset,
                                  CRSMatrixView< real64, globalIndex const > const & matrix,
                                  arrayView1d< real64 > const & rhs,
                                  LAMBDA && lambda ) const
{
  integer const component = GetComponent();
  string const & functionName = getReference< string >( viewKeyStruct::functionNameString );
  FunctionManager & functionManager = FunctionManager::Instance();

  array1d< globalIndex > dofArray( targetSet.size() );
  arrayView1d< globalIndex > const & dof = dofArray.toView();

  array1d< real64 > rhsContributionArray( targetSet.size() );
  arrayView1d< real64 > const & rhsContribution = rhsContributionArray.toView();

  real64 sizeScalingFactor = 0.0;
  if( m_normalizeBySetSize )
  {
    // note: this assumes that the ghost elements have been filtered out

    // recompute the set size here to make sure that topology changes are accounted for
    integer const localSetSize = targetSet.size();
    integer globalSetSize = 0;

    // synchronize
    MpiWrapper::allReduce( &localSetSize, &globalSetSize, 1, MPI_SUM, MPI_COMM_GEOSX );

    // set the scaling factor
    sizeScalingFactor = globalSetSize >= 1 ? 1.0 / globalSetSize : 1;
  }
  else
  {
    sizeScalingFactor = 1;
  }


  if( functionName.empty() || functionManager.getGroupReference< FunctionBase >( functionName ).isFunctionOfTime() == 2 )
  {
    real64 value = m_scale * dt * sizeScalingFactor;
    if( !functionName.empty() )
    {
      FunctionBase const & function = functionManager.getGroupReference< FunctionBase >( functionName );
      value *= function.Evaluate( &time );
    }

    forAll< POLICY >( targetSet.size(),
                      [targetSet, dof, dofMap, dofRankOffset, component, matrix, rhsContribution, value, lambda] GEOSX_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.getGroupReference< 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 * sizeScalingFactor;

    forAll< POLICY >( targetSet.size(),
                      [targetSet, dof, dofMap, dofRankOffset, component, matrix, rhsContribution, results, value, lambda] GEOSX_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 ) );
    } );
  }

  FIELD_OP::template PrescribeRhsValues< POLICY >( rhs, dof, dofRankOffset, rhsContribution );
}

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();
  forAll< POLICY >( targetSet.size(), [targetSet, dofMap, matrix, component] GEOSX_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;
    }
  } );
}

}

#endif