MimeticInnerProductHelpers.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_MIMETICINNERPRODUCTS_MIMETICINNERPRODUCTHELPERS_HPP
#define GEOS_FINITEVOLUME_MIMETICINNERPRODUCTS_MIMETICINNERPRODUCTHELPERS_HPP
 
#include "mesh/utilities/ComputationalGeometry.hpp"
 
namespace geos
{
namespace mimeticInnerProduct
{
 
struct MimeticInnerProductHelpers
{
 
  GEOS_HOST_DEVICE
  static
  void makeFullTensor( real64 const (&values)[ 3 ],
                       real64 (& result)[ 3 ][ 3 ] )
  {
    LvArray::tensorOps::fill< 3, 3 >( result, 0.0 );
    result[ 0 ][ 0 ] = values[ 0 ];
    result[ 1 ][ 1 ] = values[ 1 ];
    result[ 2 ][ 2 ] = values[ 2 ];
  }
 
  GEOS_HOST_DEVICE
  static
  void computeCellToFacetVector( real64 (& cellToFacetVec)[ 3 ],
                                 real64 const (&facetCenter)[ 3 ],
                                 arraySlice1d< real64 const > const & cellCenter )
  {
    LvArray::tensorOps::copy< 3 >( cellToFacetVec, facetCenter );
    LvArray::tensorOps::subtract< 3 >( cellToFacetVec, cellCenter );
  }
 
  GEOS_HOST_DEVICE
  static
  void orientNormalOutward( real64 const (&cellToFacetVec)[ 3 ],
                            real64 (& faceNormal)[ 3 ] )
  {
    if( LvArray::tensorOps::AiBi< 3 >( cellToFacetVec, faceNormal ) < 0.0 )
    {
      LvArray::tensorOps::scale< 3 >( faceNormal, -1.0 );
    }
  }
 
  template< localIndex NF >
  GEOS_HOST_DEVICE
  static
  void computeCellToFaceGeometry( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition,
                                  ArrayOfArraysView< localIndex const > const & faceToNodes,
                                  arraySlice1d< localIndex const > const & elemToFaces,
                                  arraySlice1d< real64 const > const & elemCenter,
                                  real64 const & areaTolerance,
                                  real64 (& C)[ NF ][ 3 ],
                                  real64 (& N)[ NF ][ 3 ],
                                  real64 (& faceArea)[ NF ] )
  {
    for( localIndex ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc )
    {
      real64 faceCenter[ 3 ], faceNormal[ 3 ], cellToFaceVec[ 3 ];
 
      faceArea[ifaceLoc] =
        computationalGeometry::centroid_3DPolygon( faceToNodes[elemToFaces[ifaceLoc]],
                                                   nodePosition,
                                                   faceCenter,
                                                   faceNormal,
                                                   areaTolerance );
 
      computeCellToFacetVector( cellToFaceVec, faceCenter, elemCenter );
      orientNormalOutward( cellToFaceVec, faceNormal );
 
      for( int d = 0; d < 3; ++d )
      {
        C[ifaceLoc][d] = cellToFaceVec[d];
        N[ifaceLoc][d] = faceArea[ifaceLoc] * faceNormal[d];
      }
    }
  }
 
  template< localIndex NF >
  GEOS_HOST_DEVICE
  static
  void computeConsistencyTerm( real64 const (&C)[ NF ][ 3 ],
                               real64 const & elemVolume,
                               real64 const (&elemPerm)[ 3 ],
                               real64 (& CKCt)[ NF ][ NF ] )
  {
    // Kinv assumes diagonal K
    real64 Kinv[ 3 ][ 3 ] = {{ 0 }};
    for( int d = 0; d < 3; ++d )
    {
      Kinv[d][d] = 1.0 / elemPerm[d];
    }
 
    real64 work[ 3 ][ NF ] = {{ 0 }};
    LvArray::tensorOps::Rij_eq_AikBjk< 3, NF, 3 >( work, Kinv, C );
    LvArray::tensorOps::Rij_eq_AikBkj< NF, NF, 3 >( CKCt, C, work );
    LvArray::tensorOps::scale< NF, NF >( CKCt, 1.0 / elemVolume );
  }
 
  template< localIndex NF >
  GEOS_HOST_DEVICE
  static
  void orthonormalize( real64 (& q0)[ NF ],
                       real64 (& q1)[ NF ],
                       real64 (& q2)[ NF ],
                       real64 (& cellToFaceMat)[ NF ][ 3 ] )
  {
    // modified Gram-Schmidt algorithm
 
    // q0
    LvArray::tensorOps::scale< NF >( q0, 1.0/LvArray::tensorOps::l2Norm< NF >( q0 ) );
 
    // q1
    real64 const q0Dotq1 = LvArray::tensorOps::AiBi< NF >( q0, q1 );
    LvArray::tensorOps::scaledAdd< NF >( q1, q0, -q0Dotq1 );
    LvArray::tensorOps::scale< NF >( q1, 1.0/LvArray::tensorOps::l2Norm< NF >( q1 ) );
 
    // q2
    real64 const q0Dotq2 = LvArray::tensorOps::AiBi< NF >( q0, q2 );
    LvArray::tensorOps::scaledAdd< NF >( q2, q0, -q0Dotq2 );
    real64 const q1Dotq2 = LvArray::tensorOps::AiBi< NF >( q1, q2 );
    LvArray::tensorOps::scaledAdd< NF >( q2, q1, -q1Dotq2 );
    LvArray::tensorOps::scale< NF >( q2, 1.0/LvArray::tensorOps::l2Norm< NF >( q2 ) );
 
    for( localIndex i = 0; i < NF; ++i )
    {
      cellToFaceMat[ i ][ 0 ] = q0[ i ];
      cellToFaceMat[ i ][ 1 ] = q1[ i ];
      cellToFaceMat[ i ][ 2 ] = q2[ i ];
    }
  }
 
  template< localIndex NF >
  GEOS_HOST_DEVICE
  static void
  computeInvTPFATransWithMultiplier( real64 const (&elemPerm)[ 3 ],
                                     real64 const (&faceNormal)[ 3 ],
                                     real64 const & faceArea,
                                     real64 const & transMult,
                                     real64 const & weightToleranceInv,
                                     real64 (& cellToFaceVec)[ 3 ],
                                     real64 & tpTransInv )
  {
    real64 const c2fDistance = LvArray::tensorOps::normalize< 3 >( cellToFaceVec );
    real64 const mult = transMult;
    tpTransInv = c2fDistance / faceArea;
 
    real64 faceConormal[ 3 ] = { 0.0 };
    LvArray::tensorOps::hadamardProduct< 3 >( faceConormal, elemPerm, faceNormal );
    real64 halfWeight = LvArray::tensorOps::AiBi< 3 >( cellToFaceVec, faceConormal );
    if( halfWeight < 0.0 )
    {
      LvArray::tensorOps::hadamardProduct< 3 >( faceConormal, elemPerm, cellToFaceVec );
      halfWeight = LvArray::tensorOps::AiBi< 3 >( cellToFaceVec, faceConormal );
    }
    tpTransInv /= halfWeight;
    tpTransInv  = LvArray::math::min( tpTransInv, weightToleranceInv );
    tpTransInv *=  ( 1.0 - mult ) / mult;
  }
 
 
  template< localIndex NF >
  GEOS_HOST_DEVICE
  static void
  computeTransMatrixWithMultipliers( real64 const (&tpTransInv)[ NF ],
                                     arraySlice2d< real64 > const & transMatrix )
  {
    // the inverse of the pertubed inverse is computed using the Sherman-Morrison formula
    for( localIndex k = 0; k < NF; ++k )
    {
      real64 const mult = LvArray::math::sqrt( tpTransInv[k] );
      real64 Tmult[ NF ] = { 0.0 };
      for( localIndex i = 0; i < NF; ++i )
      {
        Tmult[i] = transMatrix[k][i] * mult;
      }
 
      real64 const invDenom = 1.0 / ( 1.0 + Tmult[k] * mult );
      for( localIndex i = 0; i < NF; ++i )
      {
        for( localIndex j = 0; j < NF; ++j )
        {
          transMatrix[i][j] -= Tmult[i]*Tmult[j]*invDenom;
        }
      }
    }
  }
 
};
 
} // namespace mimeticInnerProduct
 
} // namespace geos
 
#endif //GEOS_FINITEVOLUME_MIMETICINNERPRODUCTS_MIMETICINNERPRODUCTHELPERS_HPP