FieldSpecificationOps.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_FIELDSPECIFICATIONOPS_HPP
#define GEOSX_MANAGERS_FIELDSPECIFICATION_FIELDSPECIFICATIONOPS_HPP

#include "common/DataTypes.hpp"

namespace geosx
{

struct OpEqual
{
  template< typename T, typename U >
  GEOSX_HOST_DEVICE static inline
  void apply( T & lhs, U const & rhs )
  {
    lhs = static_cast< T >( rhs );
  }
};

struct OpAdd
{
  template< typename T, typename U >
  GEOSX_HOST_DEVICE static inline
  void apply( T & lhs, U const & rhs )
  {
    lhs += static_cast< T >( rhs );
  }
};

template< typename OP >
struct FieldSpecificationOp
{
  using OpType = OP;

  template< typename T >
  GEOSX_HOST_DEVICE
  static inline typename std::enable_if< !traits::is_tensorT< T >, void >::type
  SpecifyFieldValue( arrayView1d< T > const & field,
                     localIndex const index,
                     integer const component,
                     real64 const value )
  {
    GEOSX_UNUSED_VAR( component );
    OP::template apply( field( index ), value );
  }

  template< typename T >
  GEOSX_HOST_DEVICE
  static inline typename std::enable_if< traits::is_tensorT< T >, void >::type
  SpecifyFieldValue( arrayView1d< T > const & field,
                     localIndex const index,
                     integer const component,
                     real64 const value )
  {
    OP::template apply( field( index )[component], value );
  }

  template< typename T >
  GEOSX_HOST_DEVICE
  static inline typename std::enable_if< !traits::is_tensorT< T >, void >::type
  ReadFieldValue( arrayView1d< T const > const & field,
                  localIndex const index,
                  integer const component,
                  real64 & value )
  {
    GEOSX_UNUSED_VAR( component );
    OP::template apply( value, field( index ) );
  }

  template< typename T >
  GEOSX_HOST_DEVICE
  static inline typename std::enable_if< traits::is_tensorT< T >, void >::type
  ReadFieldValue( arrayView1d< T const > const & field,
                  localIndex const index,
                  integer const component,
                  real64 & value )
  {
    OP::template apply( value, field( index )[component] );
  }

  template< typename T, int USD >
  GEOSX_HOST_DEVICE
  static inline typename std::enable_if< !traits::is_tensorT< T >, void >::type
  SpecifyFieldValue( arrayView2d< T, USD > const & field,
                     localIndex const index,
                     integer const component,
                     real64 const value )
  {
    if( component >= 0 )
    {
      OP::template apply( field( index, component ), value );
    }
    else
    {
      for( localIndex a = 0; a < field.size( 1 ); ++a )
      {
        OP::template apply( field( index, a ), value );
      }
    }
  }

  template< typename T, int USD >
  GEOSX_HOST_DEVICE
  static inline typename std::enable_if< traits::is_tensorT< T >, void >::type
  SpecifyFieldValue( arrayView2d< T, USD > const & field,
                     localIndex const index,
                     integer const component,
                     real64 const value )
  {
    if( component >= 0 )
    {
      for( localIndex a = 0; a < field.size( 1 ); ++a )
      {
        OP::template apply( field( index, a )[component], value );
      }
    }
    else
    {
      for( localIndex a = 0; a < field.size( 1 ); ++a )
      {
        for( localIndex c = 0; c < T::SIZE; ++c )
        {
          OP::template apply( field( index, a )[c], value );
        }
      }
    }
  }

  template< typename T, int USD >
  GEOSX_HOST_DEVICE
  static inline typename std::enable_if< !traits::is_tensorT< T >, void >::type
  ReadFieldValue( arrayView2d< T const, USD > const & field,
                  localIndex const index,
                  integer const component,
                  real64 & value )
  {
    GEOSX_ASSERT( component >= 0 );
    OP::template apply( value, field( index, component ) );
  }

  template< typename T, int USD >
  GEOSX_HOST_DEVICE
  static inline typename std::enable_if< traits::is_tensorT< T >, void >::type
  ReadFieldValue( arrayView2d< T const, USD > const & field,
                  localIndex const index,
                  integer const component,
                  real64 & value )
  {
    GEOSX_UNUSED_VAR( field );
    GEOSX_UNUSED_VAR( index );
    GEOSX_UNUSED_VAR( component );
    GEOSX_UNUSED_VAR( value );
    GEOSX_ERROR( "ReadFieldValue: unsupported operation" );
  }

  template< typename T, int USD >
  GEOSX_HOST_DEVICE
  static inline typename std::enable_if< !traits::is_tensorT< T >, void >::type
  SpecifyFieldValue( arrayView3d< T, USD > const & field,
                     localIndex const index,
                     integer const component,
                     real64 const value )
  {
    if( component >= 0 )
    {
      for( localIndex a = 0; a < field.size( 1 ); ++a )
      {
        OP::template apply( field( index, a, component ), value );
      }
    }
    else
    {
      for( localIndex a = 0; a < field.size( 1 ); ++a )
      {
        for( localIndex b = 0; b < field.size( 2 ); ++b )
        {
          OP::template apply( field( index, a, b ), value );
        }
      }
    }
  }

  template< typename T, int USD >
  GEOSX_HOST_DEVICE
  static inline typename std::enable_if< traits::is_tensorT< T >, void >::type
  SpecifyFieldValue( arrayView3d< T, USD > const & field,
                     localIndex const index,
                     integer const component,
                     real64 const value )
  {
    if( component >= 0 )
    {
      for( localIndex a = 0; a < field.size( 1 ); ++a )
      {
        for( localIndex b = 0; b < field.size( 2 ); ++b )
        {
          OP::template apply( field( index, a, b )[component], value );
        }
      }
    }
    else
    {
      for( localIndex a = 0; a < field.size( 1 ); ++a )
      {
        for( localIndex b = 0; b < field.size( 2 ); ++b )
        {
          for( localIndex c = 0; c < T::size(); ++c )
          {
            OP::template apply( field( index, a, b )[c], value );
          }
        }
      }
    }
  }

  template< typename T, int USD >
  GEOSX_HOST_DEVICE
  static inline void
  ReadFieldValue( arrayView3d< T const, USD > const & field,
                  localIndex const index,
                  integer const component,
                  real64 & value )
  {
    GEOSX_UNUSED_VAR( field );
    GEOSX_UNUSED_VAR( index );
    GEOSX_UNUSED_VAR( component );
    GEOSX_UNUSED_VAR( value );
    GEOSX_ERROR( "ReadFieldValue: unsupported operation" );
  }

};

struct FieldSpecificationEqual : public FieldSpecificationOp< OpEqual >
{
  using base_type = FieldSpecificationOp< OpEqual >;
  using base_type::SpecifyFieldValue;

  template< typename LAI >
  static inline void SpecifyFieldValue( globalIndex const dof,
                                        typename LAI::ParallelMatrix & matrix,
                                        real64 & rhs,
                                        real64 const bcValue,
                                        real64 const fieldValue )
  {
    if( matrix.getLocalRowID( dof ) >= 0 )
    {
      real64 const diag = matrix.clearRow( dof, true );
      rhs = -diag * (bcValue - fieldValue);
    }
    else
    {
      rhs = 0.0;
    }
  }

  GEOSX_HOST_DEVICE
  static inline void
  SpecifyFieldValue( globalIndex const dof,
                     globalIndex const dofRankOffset,
                     CRSMatrixView< real64, globalIndex const > const & matrix,
                     real64 & rhs,
                     real64 const bcValue,
                     real64 const fieldValue )
  {
    globalIndex const localRow = dof - dofRankOffset;
    if( localRow >= 0 && localRow < matrix.numRows() )
    {
      arraySlice1d< globalIndex const > const columns = matrix.getColumns( localRow );
      arraySlice1d< real64 > const entries = matrix.getEntries( localRow );
      localIndex const numEntries = matrix.numNonZeros( localRow );

      real64 diagonal = 0;
      for( localIndex j = 0; j < numEntries; ++j )
      {
        if( columns[ j ] == dof )
        { diagonal = entries[ j ]; }
        else
        { entries[ j ] = 0; }
      }

      rhs = -diagonal * (bcValue - fieldValue);
    }
    else
    {
      rhs = 0.0;
    }
  }

  template< typename LAI >
  static inline void PrescribeRhsValues( typename LAI::ParallelVector & rhs,
                                         localIndex const num,
                                         globalIndex * const dof,
                                         real64 * const values )
  {
    for( localIndex a = 0; a < num; ++a )
    {
      if( rhs.getLocalRowID( dof[a] ) >= 0 )
      {
        rhs.set( dof[a], values[a] );
      }
    }
  }

  template< typename POLICY >
  static inline void PrescribeRhsValues( arrayView1d< real64 > const & rhs,
                                         arrayView1d< globalIndex const > const & dof,
                                         globalIndex const dofRankOffset,
                                         arrayView1d< real64 const > const & values )
  {
    GEOSX_ASSERT_EQ( dof.size(), values.size() );
    forAll< POLICY >( dof.size(), [rhs, dof, dofRankOffset, values] GEOSX_HOST_DEVICE ( localIndex const a )
    {
      globalIndex const localRow = dof[ a ] - dofRankOffset;
      if( localRow >= 0 && localRow < rhs.size() )
      { rhs[ localRow ] = values[ a ]; }
    } );
  }
};

struct FieldSpecificationAdd : public FieldSpecificationOp< OpAdd >
{
  using base_type = FieldSpecificationOp< OpAdd >;
  using base_type::SpecifyFieldValue;

  GEOSX_HOST_DEVICE
  static inline void
  SpecifyFieldValue( globalIndex const dof,
                     globalIndex const dofRankOffset,
                     CRSMatrixView< real64, globalIndex const > const & matrix,
                     real64 & rhs,
                     real64 const bcValue,
                     real64 const fieldValue )
  {
    GEOSX_UNUSED_VAR( dof );
    GEOSX_UNUSED_VAR( dofRankOffset );
    GEOSX_UNUSED_VAR( matrix );
    GEOSX_UNUSED_VAR( fieldValue );
    rhs += bcValue;
  }

  template< typename LAI >
  static inline void PrescribeRhsValues( typename LAI::ParallelVector & rhs,
                                         localIndex const num,
                                         globalIndex * const dof,
                                         real64 * const values )
  {
    rhs.add( dof, values, num );
  }

  template< typename POLICY >
  static inline void PrescribeRhsValues( arrayView1d< real64 > const & rhs,
                                         arrayView1d< globalIndex const > const & dof,
                                         globalIndex const dofRankOffset,
                                         arrayView1d< real64 const > const & values )
  {
    GEOSX_ASSERT_EQ( dof.size(), values.size() );
    forAll< POLICY >( dof.size(), [rhs, dof, dofRankOffset, values] GEOSX_HOST_DEVICE ( localIndex const a )
    {
      globalIndex const localRow = dof[ a ] - dofRankOffset;
      if( localRow >= 0 && localRow < rhs.size() )
      { rhs[ localRow ] += values[ a ]; }
    } );
  }

};

} //namespace geosx

#endif //GEOSX_MANAGERS_FIELDSPECIFICATION_FIELDSPECIFICATIONOPS_HPP