BoundaryStencil.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 TotalEnergies
 * Copyright (c) 2019-     GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */
 
#ifndef GEOS_FINITEVOLUME_BOUNDARYSTENCIL_HPP_
#define GEOS_FINITEVOLUME_BOUNDARYSTENCIL_HPP_
 
#include "StencilBase.hpp"
 
namespace geos
{
 
class BoundaryStencilWrapper : public StencilWrapperBase< TwoPointStencilTraits >
{
public:
 
  template< typename VIEWTYPE >
  using CoefficientAccessor = ElementRegionManager::ElementViewConst< VIEWTYPE >;
 
  BoundaryStencilWrapper( IndexContainerType const & elementRegionIndices,
                          IndexContainerType const & elementSubRegionIndices,
                          IndexContainerType const & elementIndices,
                          WeightContainerType const & weights,
                          arrayView2d< real64, nodes::REFERENCE_POSITION_USD > const & faceNormal,
                          arrayView2d< real64, nodes::REFERENCE_POSITION_USD > const & cellToFaceVec,
                          arrayView1d< real64 > const & weightMultiplier );
 
  GEOS_HOST_DEVICE
  void computeWeights( localIndex const iconn,
                       CoefficientAccessor< arrayView3d< real64 const > > const & coefficient,
                       real64 & weight,
                       real64 ( &dWeight_dCoef )[3] ) const;
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  localIndex size() const
  { return m_elementRegionIndices.size( 0 ); }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr localIndex stencilSize( localIndex const index ) const
  {
    GEOS_UNUSED_VAR( index );
    return maxStencilSize;
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  constexpr localIndex numPointsInFlux( localIndex const index ) const
  {
    GEOS_UNUSED_VAR( index );
    return maxNumPointsInFlux;
  }
 
private:
 
  arrayView2d< real64, nodes::REFERENCE_POSITION_USD > m_faceNormal;
  arrayView2d< real64, nodes::REFERENCE_POSITION_USD > m_cellToFaceVec;
  arrayView1d< real64 > m_weightMultiplier;
};
 
class BoundaryStencil final : public StencilBase< TwoPointStencilTraits, BoundaryStencil >
{
public:
 
  BoundaryStencil();
 
  struct Order
  {
    static constexpr localIndex ELEM = 0; 
    static constexpr localIndex FACE = 1; 
  };
 
  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] );
 
  virtual localIndex size() const override
  {
    return m_elementRegionIndices.size( 0 );
  }
 
  constexpr localIndex stencilSize( localIndex const index ) const
  {
    GEOS_UNUSED_VAR( index );
    return maxStencilSize;
  }
 
  using KernelWrapper = BoundaryStencilWrapper;
 
  KernelWrapper createKernelWrapper() const;
 
private:
 
  array2d< real64, nodes::REFERENCE_POSITION_PERM > m_faceNormal;
  array2d< real64, nodes::REFERENCE_POSITION_PERM > m_cellToFaceVec;
  array1d< real64 > m_weightMultiplier;
 
};
 
GEOS_HOST_DEVICE
inline void
BoundaryStencilWrapper::
  computeWeights( localIndex const iconn,
                  CoefficientAccessor< arrayView3d< real64 const > > const & coefficient,
                  real64 & weight,
                  real64 (& dWeight_dCoef)[3] ) const
{
  localIndex const er  = m_elementRegionIndices[iconn][BoundaryStencil::Order::ELEM];
  localIndex const esr = m_elementSubRegionIndices[iconn][BoundaryStencil::Order::ELEM];
  localIndex const ei  = m_elementIndices[iconn][BoundaryStencil::Order::ELEM];
 
  // Preload data onto the stack
  real64 faceNormal[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( m_faceNormal[iconn] );
  real64 const cellToFace[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( m_cellToFaceVec[iconn] );
  arraySlice1d< real64 const > const coef = coefficient[er][esr][ei][0];
 
  // Compute the face-cell transmissibility
  auto const compute = [&]
  {
    real64 faceConormal[3];
    LvArray::tensorOps::hadamardProduct< 3 >( faceConormal, coef, faceNormal );
    weight = LvArray::tensorOps::AiBi< 3 >( cellToFace, faceConormal );
    LvArray::tensorOps::hadamardProduct< 3 >( dWeight_dCoef, cellToFace, faceNormal );
  };
 
  compute();
  if( weight < 0.0 )
  {
    // Correct negative weight issue arising from non-K-orthogonal grids
    // by replacing face normal with c2f vector
    LvArray::tensorOps::copy< 3 >( faceNormal, cellToFace );
    compute();
  }
 
  // Scale by face/distance and trans multiplier
  real64 const mult = m_weights[iconn][BoundaryStencil::Order::ELEM] * m_weightMultiplier[iconn];
  weight *= mult;
  LvArray::tensorOps::scale< 3 >( dWeight_dCoef, mult );
}
 
} // namespace geos
 
#endif // GEOS_FINITEVOLUME_BOUNDARYSTENCIL_HPP_