FieldSpecificationOps.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_COMMON_FIELDSPECIFICATIONOPS_HPP
#define GEOS_COMMON_FIELDSPECIFICATIONOPS_HPP
 
#include "codingUtilities/traits.hpp"
#include "common/DataTypes.hpp"
#include "common/GEOS_RAJA_Interface.hpp"
 
namespace geos
{
 
struct OpEqual
{
  template< typename T, typename U >
  GEOS_HOST_DEVICE static inline
  void apply( T & lhs, U const & rhs )
  {
    lhs = static_cast< T >( rhs );
  }
};
 
struct OpAdd
{
  template< typename T, typename U >
  GEOS_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 >
  GEOS_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 )
  {
    GEOS_UNUSED_VAR( component );
    OP::template apply( field( index ), value );
  }
 
  template< typename T >
  GEOS_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 >
  GEOS_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 )
  {
    GEOS_UNUSED_VAR( component );
    OP::template apply( value, field( index ) );
  }
 
  template< typename T >
  GEOS_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 >
  GEOS_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 >
  GEOS_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 >
  GEOS_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 )
  {
    GEOS_ASSERT( component >= 0 );
    OP::template apply( value, field( index, component ) );
  }
 
  template< typename T, int USD >
  GEOS_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 )
  {
    GEOS_UNUSED_VAR( field );
    GEOS_UNUSED_VAR( index );
    GEOS_UNUSED_VAR( component );
    GEOS_UNUSED_VAR( value );
    GEOS_ERROR( "readFieldValue: unsupported operation" );
  }
 
  template< typename T, int USD >
  GEOS_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 >
  GEOS_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 >
  GEOS_HOST_DEVICE
  static inline void
  readFieldValue( arrayView3d< T const, USD > const & field,
                  localIndex const index,
                  integer const component,
                  real64 & value )
  {
    GEOS_UNUSED_VAR( field );
    GEOS_UNUSED_VAR( index );
    GEOS_UNUSED_VAR( component );
    GEOS_UNUSED_VAR( value );
    GEOS_ERROR( "readFieldValue: unsupported operation" );
  }
 
  template< typename T, int USD >
  GEOS_HOST_DEVICE
  static inline typename std::enable_if< !traits::is_tensorT< T >, void >::type
  SpecifyFieldValue( arrayView4d< 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 < field.size( 3 ); ++c )
          {
            OP::template apply( field( index, a, b, c ), value );
          }
        }
      }
    }
  }
 
  template< typename T, int USD >
  GEOS_HOST_DEVICE
  static inline typename std::enable_if< traits::is_tensorT< T >, void >::type
  SpecifyFieldValue( arrayView4d< 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 )
        {
          for( localIndex c = 0; c < field.size( 3 ); ++c )
          {
            OP::template apply( field( index, a, b, c )[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 < field.size( 3 ); ++c )
          {
            for( localIndex d = 0; d < T::size(); ++d )
            {
              OP::template apply( field( index, a, b, c )[d], value );
            }
          }
        }
      }
    }
  }
 
  template< typename T, int USD >
  GEOS_HOST_DEVICE
  static inline void
  readFieldValue( arrayView4d< T const, USD > const & field,
                  localIndex const index,
                  integer const component,
                  real64 & value )
  {
    GEOS_UNUSED_VAR( field );
    GEOS_UNUSED_VAR( index );
    GEOS_UNUSED_VAR( component );
    GEOS_UNUSED_VAR( value );
    GEOS_ERROR( "readFieldValue: unsupported operation" );
  }
 
};
 
struct FieldSpecificationEqual : public FieldSpecificationOp< OpEqual >
{
  using base_type = FieldSpecificationOp< OpEqual >;
  using base_type::SpecifyFieldValue;
 
  GEOS_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;
      real64 const minDiagonal = 1e-15;
      for( localIndex j = 0; j < numEntries; ++j )
      {
        if( columns[ j ] == dof )
        {
          // check that the entry is large enough to enforce the boundary condition
          if( entries[j] >= 0 && entries[j] < minDiagonal )
          {
            entries[ j ] = minDiagonal;
          }
          else if( entries[j] < 0 && entries[j] > -minDiagonal )
          {
            entries[ j ] = -minDiagonal;
          }
          diagonal = entries[j];
        }
        else
        { entries[ j ] = 0; }
      }
 
      rhs = -diagonal * (bcValue - fieldValue);
    }
    else
    {
      rhs = 0.0;
    }
  }
 
  template< typename POLICY >
  static inline void prescribeRhsValues( arrayView1d< real64 > const & rhs,
                                         arrayView1d< globalIndex const > const & dof,
                                         globalIndex const dofRankOffset,
                                         arrayView1d< real64 const > const & values )
  {
    GEOS_ASSERT_EQ( dof.size(), values.size() );
    forAll< POLICY >( dof.size(), [rhs, dof, dofRankOffset, values] GEOS_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;
 
  GEOS_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 )
  {
    GEOS_UNUSED_VAR( dof );
    GEOS_UNUSED_VAR( dofRankOffset );
    GEOS_UNUSED_VAR( matrix );
    GEOS_UNUSED_VAR( fieldValue );
    rhs += bcValue;
  }
 
  template< typename POLICY >
  static inline void prescribeRhsValues( arrayView1d< real64 > const & rhs,
                                         arrayView1d< globalIndex const > const & dof,
                                         globalIndex const dofRankOffset,
                                         arrayView1d< real64 const > const & values )
  {
    GEOS_ASSERT_EQ( dof.size(), values.size() );
    forAll< POLICY >( dof.size(), [rhs, dof, dofRankOffset, values] GEOS_HOST_DEVICE ( localIndex const a )
    {
      globalIndex const localRow = dof[ a ] - dofRankOffset;
      if( localRow >= 0 && localRow < rhs.size() )
      { rhs[ localRow ] += values[ a ]; }
    } );
  }
 
};
 
} //namespace geos
 
#endif //GEOS_COMMON_FIELDSPECIFICATIONOPS_HPP