FaceElementToCellStencil.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_FINITEVOLUME_FACEELEMENTTOCELLSTENCIL_HPP_
#define GEOS_FINITEVOLUME_FACEELEMENTTOCELLSTENCIL_HPP_
 
#include "StencilBase.hpp"
 
namespace geos
{
 
class FaceElementToCellStencilWrapper : public StencilWrapperBase< TwoPointStencilTraits >
{
public:
 
  template< typename VIEWTYPE >
  using CoefficientAccessor = ElementRegionManager::ElementViewConst< VIEWTYPE >;
 
  FaceElementToCellStencilWrapper( IndexContainerType const & elementRegionIndices,
                                   IndexContainerType const & elementSubRegionIndices,
                                   IndexContainerType const & elementIndices,
                                   WeightContainerType const & weights,
                                   arrayView2d< real64 > const & faceNormal,
                                   arrayView2d< real64 > const & cellToFaceVec,
                                   arrayView1d< real64 > const & transMultiplier );
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  localIndex size() const
  {
    return m_elementRegionIndices.size( 0 );
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr localIndex stencilSize( localIndex index ) const
  {
    GEOS_UNUSED_VAR( index );
    return maxStencilSize;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr localIndex numPointsInFlux( localIndex index ) const
  {
    GEOS_UNUSED_VAR( index );
    return maxNumPointsInFlux;
  }
 
  GEOS_HOST_DEVICE
  void computeWeights( localIndex iconn,
                       CoefficientAccessor< arrayView3d< real64 const > > const & coefficient,
                       CoefficientAccessor< arrayView3d< real64 const > > const & dCoeff_dVar,
                       real64 ( &weight )[1][2],
                       real64 ( &dWeight_dVar )[1][2] ) const;
 
  GEOS_HOST_DEVICE
  void computeWeights( localIndex iconn,
                       real64 ( &weight )[1][2],
                       real64 ( &dWeight_dVar )[1][2] ) const;
 
  GEOS_HOST_DEVICE
  void computeWeights( localIndex iconn,
                       CoefficientAccessor< arrayView3d< real64 const > > const & coefficient,
                       CoefficientAccessor< arrayView3d< real64 const > > const & dCoeff_dVar1,
                       CoefficientAccessor< arrayView3d< real64 const > > const & dCoeff_dVar2,
                       real64 ( &weight )[1][2],
                       real64 ( &dWeight_dVar1 )[1][2],
                       real64 ( &dWeight_dVar2 )[1][2] ) const;
 
  GEOS_HOST_DEVICE
  void computeStabilizationWeights( localIndex iconn,
                                    real64 ( & stabilizationWeight )[1][2] ) const
  { GEOS_UNUSED_VAR( iconn, stabilizationWeight ); }
 
  GEOS_HOST_DEVICE
  void removeHydraulicApertureContribution( localIndex const iconn, ElementRegionManager::ElementViewConst< arrayView1d< real64 const > > hydraulicAperture ) const;
 
  GEOS_HOST_DEVICE
  void addHydraulicApertureContribution( localIndex const iconn, ElementRegionManager::ElementViewConst< arrayView1d< real64 const > > hydraulicAperture ) const;
 
private:
 
  arrayView2d< real64 > m_faceNormal;
 
  arrayView2d< real64 > m_cellToFaceVec;
 
  arrayView1d< real64 > m_transMultiplier;
};
 
class FaceElementToCellStencil final : public StencilBase< TwoPointStencilTraits, FaceElementToCellStencil >
{
public:
 
  FaceElementToCellStencil();
 
  virtual void move( LvArray::MemorySpace const space ) override;
 
  virtual void add( localIndex const numPts,
                    localIndex const * const elementRegionIndices,
                    localIndex const * const elementSubRegionIndices,
                    localIndex const * const elementIndices,
                    real64 const * const weights,
                    localIndex const connectorIndex ) override;
 
  void addVectors( real64 const & transMultiplier,
                   real64 const (&faceNormal)[3],
                   real64 const (&cellToFaceVec)[3] );
 
  using KernelWrapper = FaceElementToCellStencilWrapper;
 
  KernelWrapper createKernelWrapper() const;
 
  virtual localIndex size() const override
  { return m_elementRegionIndices.size( 0 ); }
 
  virtual void reserve( localIndex const size ) override;
 
  constexpr localIndex stencilSize( localIndex index ) const
  {
    GEOS_UNUSED_VAR( index );
    return maxStencilSize;
  }
 
private:
 
  array2d< real64 > m_faceNormal;
  array2d< real64 > m_cellToFaceVec;
  array1d< real64 > m_transMultiplier;
};
 
GEOS_HOST_DEVICE
inline void FaceElementToCellStencilWrapper::
  computeWeights( localIndex const iconn,
                  CoefficientAccessor< arrayView3d< real64 const > > const & coefficient,
                  CoefficientAccessor< arrayView3d< real64 const > > const & dCoeff_dVar,
                  real64 ( & weight )[1][2],
                  real64 ( & dWeight_dVar )[1][2] ) const
{
  localIndex const er0  =  m_elementRegionIndices[iconn][0];
  localIndex const esr0 =  m_elementSubRegionIndices[iconn][0];
  localIndex const ei0  =  m_elementIndices[iconn][0];
 
  localIndex const er1  =  m_elementRegionIndices[iconn][1];
  localIndex const esr1 =  m_elementSubRegionIndices[iconn][1];
  localIndex const ei1  =  m_elementIndices[iconn][1];
 
  real64 faceConormal[3];
 
  // Will change when implementing collocation points.
  LvArray::tensorOps::hadamardProduct< 3 >( faceConormal, coefficient[er0][esr0][ei0][0], m_faceNormal[iconn] );
  real64 t0 = m_weights[iconn][0] * LvArray::tensorOps::AiBi< 3 >( m_cellToFaceVec[iconn], faceConormal );
  if( t0 < 0.0 )
  {
    LvArray::tensorOps::hadamardProduct< 3 >( faceConormal,
                                              coefficient[er0][esr0][ei0][0],
                                              m_cellToFaceVec[iconn] );
    t0 = m_weights[iconn][0] * LvArray::tensorOps::AiBi< 3 >( m_cellToFaceVec[iconn], faceConormal );
  }
 
  // We consider the 3rd component of the permeability which is the normal one.
  real64 const t1 = m_weights[iconn][1] * coefficient[er1][esr1][ei1][0][2];
 
  real64 const sumOfTrans = t0+t1;
  real64 const value = m_transMultiplier[iconn]*t0*t1/sumOfTrans;
 
  weight[0][0] = value;
  weight[0][1] = -value;
 
  // We consider the 3rd component of the permeability which is the normal one.
  real64 const dt0 = m_weights[iconn][0] * dCoeff_dVar[er0][esr0][ei0][0][0];
  real64 const dt1 = m_weights[iconn][1] * dCoeff_dVar[er1][esr1][ei1][0][2];
 
  dWeight_dVar[0][0] = m_transMultiplier[iconn] * ( dt0 * t1 * sumOfTrans - dt0 * t0 * t1 ) / ( sumOfTrans * sumOfTrans );
  dWeight_dVar[0][1] = m_transMultiplier[iconn] * ( t0 * dt1 * sumOfTrans - dt1 * t0 * t1 ) / ( sumOfTrans * sumOfTrans );
}
 
GEOS_HOST_DEVICE
inline void
FaceElementToCellStencilWrapper
  ::computeWeights( localIndex iconn,
                    real64 ( & weight )[1][2],
                    real64 ( & dWeight_dVar )[1][2] ) const
{
  // Will change when implementing collocation points.
  real64 const t0 = m_weights[iconn][0] * LvArray::tensorOps::AiBi< 3 >( m_cellToFaceVec[iconn], m_faceNormal[iconn] );
  real64 const t1 = m_weights[iconn][1];
 
  real64 const sumOfTrans = t0+t1;
  real64 const value = m_transMultiplier[iconn]*t0*t1/sumOfTrans;
 
  weight[0][0] = value;
  weight[0][1] = -value;
 
  dWeight_dVar[0][0] = 0.0;
  dWeight_dVar[0][1] = 0.0;
}
 
GEOS_HOST_DEVICE
inline void
FaceElementToCellStencilWrapper::
  computeWeights( localIndex const iconn,
                  CoefficientAccessor< arrayView3d< real64 const > > const & coefficient,
                  CoefficientAccessor< arrayView3d< real64 const > > const & dCoeff_dVar1,
                  CoefficientAccessor< arrayView3d< real64 const > > const & dCoeff_dVar2,
                  real64 (& weight)[1][2],
                  real64 (& dWeight_dVar1 )[1][2],
                  real64 (& dWeight_dVar2 )[1][2] ) const
{
  localIndex const er0  =  m_elementRegionIndices[iconn][0];
  localIndex const esr0 =  m_elementSubRegionIndices[iconn][0];
  localIndex const ei0  =  m_elementIndices[iconn][0];
 
  localIndex const er1  =  m_elementRegionIndices[iconn][1];
  localIndex const esr1 =  m_elementSubRegionIndices[iconn][1];
  localIndex const ei1  =  m_elementIndices[iconn][1];
 
  real64 faceConormal[3];
 
  // Will change when implementing collocation points.
  LvArray::tensorOps::hadamardProduct< 3 >( faceConormal, coefficient[er0][esr0][ei0][0], m_faceNormal[iconn] );
  real64 t0 = m_weights[iconn][0] * LvArray::tensorOps::AiBi< 3 >( m_cellToFaceVec[iconn], faceConormal );
  if( t0 < 0.0 )
  {
    LvArray::tensorOps::hadamardProduct< 3 >( faceConormal,
                                              coefficient[er0][esr0][ei0][0],
                                              m_cellToFaceVec[iconn] );
    t0 = m_weights[iconn][0] * LvArray::tensorOps::AiBi< 3 >( m_cellToFaceVec[iconn], faceConormal );
  }
  // We consider the 3rd component of the permeability which is the normal one.
  real64 const t1 = m_weights[iconn][1] * coefficient[er1][esr1][ei1][0][2];
 
  real64 const sumOfTrans = t0+t1;
  real64 const value = m_transMultiplier[iconn]*t0*t1/sumOfTrans;
 
  weight[0][0] = value;
  weight[0][1] = -value;
 
  // We consider the 3rd component of the permeability which is the normal one.
  real64 const dt0_dVar1 = m_weights[iconn][0] * dCoeff_dVar1[er0][esr0][ei0][0][0];
  real64 const dt1_dVar1 = m_weights[iconn][1] * dCoeff_dVar1[er1][esr1][ei1][0][2];
  real64 const dt0_dVar2 = m_weights[iconn][0] * dCoeff_dVar2[er0][esr0][ei0][0][0];
  real64 const dt1_dVar2 = m_weights[iconn][1] * dCoeff_dVar2[er1][esr1][ei1][0][2];
 
  dWeight_dVar1[0][0] = m_transMultiplier[iconn] * ( dt0_dVar1 * t1 * sumOfTrans - dt0_dVar1 * t0 * t1 ) / ( sumOfTrans * sumOfTrans );
  dWeight_dVar1[0][1] = m_transMultiplier[iconn] * ( t0 * dt1_dVar1 * sumOfTrans - dt1_dVar1 * t0 * t1 ) / ( sumOfTrans * sumOfTrans );
 
  dWeight_dVar2[0][0] = m_transMultiplier[iconn] * ( dt0_dVar2 * t1 * sumOfTrans - dt0_dVar2 * t0 * t1 ) / ( sumOfTrans * sumOfTrans );
  dWeight_dVar2[0][1] = m_transMultiplier[iconn] * ( t0 * dt1_dVar2 * sumOfTrans - dt1_dVar2 * t0 * t1 ) / ( sumOfTrans * sumOfTrans );
}
 
GEOS_HOST_DEVICE
inline void
FaceElementToCellStencilWrapper::
  removeHydraulicApertureContribution( localIndex const iconn, ElementRegionManager::ElementViewConst< arrayView1d< real64 const > > hydraulicAperture ) const
{
  // only the fracture side is modified, k=1
  localIndex constexpr k = 1;
  localIndex const er  =  m_elementRegionIndices[iconn][k];
  localIndex const esr =  m_elementSubRegionIndices[iconn][k];
  localIndex const ei  =  m_elementIndices[iconn][k];
 
  m_weights[iconn][k] = m_weights[iconn][k] * hydraulicAperture[er][esr][ei];
}
 
GEOS_HOST_DEVICE
inline void
FaceElementToCellStencilWrapper::
  addHydraulicApertureContribution( localIndex const iconn, ElementRegionManager::ElementViewConst< arrayView1d< real64 const > > hydraulicAperture ) const
{
  // only the fracture side is modified, k=1
  localIndex constexpr k = 1;
  localIndex const er  =  m_elementRegionIndices[iconn][k];
  localIndex const esr =  m_elementSubRegionIndices[iconn][k];
  localIndex const ei  =  m_elementIndices[iconn][k];
 
  m_weights[iconn][k] = m_weights[iconn][k] / hydraulicAperture[er][esr][ei];
}
 
} /* namespace geos */
 
#endif /* GEOS_FINITEVOLUME_FACEELEMENTTOCELLSTENCIL_HPP_ */