KernelBase.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_FINITEELEMENT_KERNELBASE_HPP_
#define GEOS_FINITEELEMENT_KERNELBASE_HPP_
 
#include "common/DataTypes.hpp"
#include "common/TimingMacros.hpp"
#include "constitutive/ConstitutivePassThru.hpp"
#include "finiteElement/FiniteElementDispatch.hpp"
#include "mesh/MeshLevel.hpp"
#include "common/GEOS_RAJA_Interface.hpp"
#include "common/TypeDispatch.hpp"
 
#include <utility>
 
#ifndef SELECTED_FE_TYPES
#define SELECTED_FE_TYPES BASE_FE_TYPES
#endif
 
namespace geos
{
 
namespace finiteElement
{
 
template< typename SUBREGION_TYPE,
          typename CONSTITUTIVE_TYPE,
          typename FE_TYPE,
          int NUM_DOF_PER_TEST_SP,
          int NUM_DOF_PER_TRIAL_SP >
class KernelBase
{
public:
  static constexpr int maxNumTestSupportPointsPerElem  = FE_TYPE::maxSupportPoints;
 
  static constexpr int maxNumTrialSupportPointsPerElem = FE_TYPE::maxSupportPoints;
 
  static constexpr int numDofPerTestSupportPoint    = NUM_DOF_PER_TEST_SP;
 
  static constexpr int numDofPerTrialSupportPoint   = NUM_DOF_PER_TRIAL_SP;
 
  static constexpr int numQuadraturePointsPerElem = FE_TYPE::numQuadraturePoints;
 
  KernelBase( SUBREGION_TYPE const & elementSubRegion,
              FE_TYPE const & finiteElementSpace,
              CONSTITUTIVE_TYPE & inputConstitutiveType ):
    m_elemsToNodes( elementSubRegion.nodeList().toViewConst() ),
    m_elemGhostRank( elementSubRegion.ghostRank() ),
    m_constitutiveUpdate( inputConstitutiveType.createKernelUpdates() ),
    m_finiteElementSpace( finiteElementSpace )
  {}
 
  struct StackVariables
  {};
 
  GEOS_HOST_DEVICE
  inline
  void setup( localIndex const k,
              StackVariables & stack ) const
  {
    GEOS_UNUSED_VAR( k );
    GEOS_UNUSED_VAR( stack );
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  void quadraturePointKernel( localIndex const k,
                              localIndex const q,
                              StackVariables & stack ) const
  {
    GEOS_UNUSED_VAR( k );
    GEOS_UNUSED_VAR( q );
    GEOS_UNUSED_VAR( stack );
  }
 
  GEOS_HOST_DEVICE
  inline
  real64 complete( localIndex const k,
                   StackVariables & stack ) const
  {
    GEOS_UNUSED_VAR( k );
    GEOS_UNUSED_VAR( stack );
    return 0;
  }
 
 
  //START_kernelLauncher
  template< typename POLICY,
            typename KERNEL_TYPE >
  static
  real64
  kernelLaunch( localIndex const numElems,
                KERNEL_TYPE const & kernelComponent )
  {
    GEOS_MARK_FUNCTION;
 
    // Define a RAJA reduction variable to get the maximum residual contribution.
    RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxResidual( 0 );
 
    forAll< POLICY >( numElems,
                      [=] GEOS_HOST_DEVICE ( localIndex const k )
    {
      typename KERNEL_TYPE::StackVariables stack;
 
      kernelComponent.setup( k, stack );
      // #pragma unroll
      for( integer q=0; q<numQuadraturePointsPerElem; ++q )
      {
        kernelComponent.quadraturePointKernel( k, q, stack );
      }
      maxResidual.max( kernelComponent.complete( k, stack ) );
    } );
    return maxResidual.get();
  }
  //END_kernelLauncher
 
protected:
  using KernelWrapper = decltype( std::declval< CONSTITUTIVE_TYPE const & >().createKernelUpdates() );
 
  traits::ViewTypeConst< typename SUBREGION_TYPE::NodeMapType::base_type > const m_elemsToNodes;
 
  arrayView1d< integer const > const m_elemGhostRank;
 
  KernelWrapper const m_constitutiveUpdate;
 
  FE_TYPE const m_finiteElementSpace;
};
 
template< template< typename SUBREGION_TYPE,
                    typename CONSTITUTIVE_TYPE,
                    typename FE_TYPE > class KERNEL_TYPE,
          typename ... ARGS >
class KernelFactory
{
public:
 
  KernelFactory( ARGS ... args ):
    m_args( args ... )
  {}
 
  template< typename SUBREGION_TYPE, typename CONSTITUTIVE_TYPE, typename FE_TYPE >
  KERNEL_TYPE< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE > createKernel(
    NodeManager & nodeManager,
    EdgeManager const & edgeManager,
    FaceManager const & faceManager,
    localIndex const targetRegionIndex,
    SUBREGION_TYPE const & elementSubRegion,
    FE_TYPE const & finiteElementSpace,
    CONSTITUTIVE_TYPE & inputConstitutiveType )
  {
    camp::tuple< NodeManager &,
                 EdgeManager const &,
                 FaceManager const &,
                 localIndex const,
                 SUBREGION_TYPE const &,
                 FE_TYPE const &,
                 CONSTITUTIVE_TYPE & > standardArgs { nodeManager,
                                                      edgeManager,
                                                      faceManager,
                                                      targetRegionIndex,
                                                      elementSubRegion,
                                                      finiteElementSpace,
                                                      inputConstitutiveType };
 
    auto allArgs = camp::tuple_cat_pair( standardArgs, m_args );
    return camp::make_from_tuple< KERNEL_TYPE< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE > >( allArgs );
  }
 
private:
  camp::tuple< ARGS ... > m_args;
};
 
 
//*****************************************************************************
//*****************************************************************************
//*****************************************************************************
 
//START_regionBasedKernelApplication
template< typename POLICY,
          typename CONSTITUTIVE_BASE,
          typename SUBREGION_TYPE,
          typename KERNEL_FACTORY >
static
real64 regionBasedKernelApplication( MeshLevel & mesh,
                                     string_array const & targetRegions,
                                     string const & finiteElementName,
                                     string const & constitutiveStringName,
                                     KERNEL_FACTORY & kernelFactory )
{
  GEOS_MARK_FUNCTION;
  // save the maximum residual contribution for scaling residuals for convergence criteria.
  real64 maxResidualContribution = 0;
 
  NodeManager & nodeManager = mesh.getNodeManager();
  EdgeManager & edgeManager = mesh.getEdgeManager();
  FaceManager & faceManager = mesh.getFaceManager();
  ElementRegionManager & elementRegionManager = mesh.getElemManager();
 
  // Loop over all sub-regions in regions of type SUBREGION_TYPE, that are listed in the targetRegions array.
  elementRegionManager.forElementSubRegions< SUBREGION_TYPE >( targetRegions,
                                                               [&constitutiveStringName,
                                                                &maxResidualContribution,
                                                                &nodeManager,
                                                                &edgeManager,
                                                                &faceManager,
                                                                &kernelFactory,
                                                                &finiteElementName]
                                                                 ( localIndex const targetRegionIndex, auto & elementSubRegion )
  {
    localIndex const numElems = elementSubRegion.size();
 
    // Get the constitutive model...and allocate a null constitutive model if required.
 
    constitutive::ConstitutiveBase * constitutiveRelation = nullptr;
    constitutive::NullModel * nullConstitutiveModel = nullptr;
    if( elementSubRegion.template hasWrapper< string >( constitutiveStringName ) )
    {
      string const & constitutiveName = elementSubRegion.template getReference< string >( constitutiveStringName );
      constitutiveRelation = &elementSubRegion.getConstitutiveModel( constitutiveName );
    }
    else
    {
      nullConstitutiveModel = &elementSubRegion.template registerGroup< constitutive::NullModel >( "nullModelGroup" );
      constitutiveRelation = nullConstitutiveModel;
    }
 
    // Call the constitutive dispatch which converts the type of constitutive model into a compile time constant.
    constitutive::ConstitutivePassThru< CONSTITUTIVE_BASE >::execute( *constitutiveRelation,
                                                                      [&maxResidualContribution,
                                                                       &nodeManager,
                                                                       &edgeManager,
                                                                       &faceManager,
                                                                       targetRegionIndex,
                                                                       &kernelFactory,
                                                                       &elementSubRegion,
                                                                       &finiteElementName,
                                                                       numElems]
                                                                        ( auto & castedConstitutiveRelation )
    {
      FiniteElementBase &
      subRegionFE = elementSubRegion.template getReference< FiniteElementBase >( finiteElementName );
 
      finiteElement::FiniteElementDispatchHandler< SELECTED_FE_TYPES >::dispatch3D( subRegionFE,
                                                                                    [&maxResidualContribution,
                                                                                     &nodeManager,
                                                                                     &edgeManager,
                                                                                     &faceManager,
                                                                                     targetRegionIndex,
                                                                                     &kernelFactory,
                                                                                     &elementSubRegion,
                                                                                     numElems,
                                                                                     &castedConstitutiveRelation] ( auto const finiteElement )
      {
        auto kernel = kernelFactory.createKernel( nodeManager,
                                                  edgeManager,
                                                  faceManager,
                                                  targetRegionIndex,
                                                  elementSubRegion,
                                                  finiteElement,
                                                  castedConstitutiveRelation );
 
        using KERNEL_TYPE = decltype( kernel );
 
        // Call the kernelLaunch function, and store the maximum contribution to the residual.
        maxResidualContribution =
          std::max( maxResidualContribution,
                    KERNEL_TYPE::template kernelLaunch< POLICY, KERNEL_TYPE >( numElems, kernel ) );
      } );
    } );
 
    // Remove the null constitutive model (not required, but cleaner)
    if( nullConstitutiveModel )
    {
      elementSubRegion.deregisterGroup( "nullModelGroup" );
    }
 
  } );
 
  return maxResidualContribution;
}
 
template< typename POLICY,
          typename DISPATCH_TYPE_LIST,
          typename KERNEL_FACTORY >
static
real64 regionBasedKernelApplication( MeshLevel & mesh,
                                     string_array const & targetRegions,
                                     string const & finiteElementName,
                                     string const & constitutiveStringName,
                                     KERNEL_FACTORY & kernelFactory )
{
  GEOS_MARK_FUNCTION;
  // save the maximum residual contribution for scaling residuals for convergence criteria.
  real64 maxResidualContribution = 0;
 
  NodeManager & nodeManager = mesh.getNodeManager();
  EdgeManager & edgeManager = mesh.getEdgeManager();
  FaceManager & faceManager = mesh.getFaceManager();
  ElementRegionManager & elementRegionManager = mesh.getElemManager();
 
  // Currently only CellElementSubRegion is supported.
  using SUBREGION_TYPES = types::TypeList< CellElementSubRegion >;
 
  // Loop over all sub-regions in regions of type SUBREGION_TYPE, that are listed in the targetRegions array.
  elementRegionManager.forElementSubRegions( SUBREGION_TYPES{}, targetRegions,
                                             [&constitutiveStringName,
                                              &maxResidualContribution,
                                              &nodeManager,
                                              &edgeManager,
                                              &faceManager,
                                              &kernelFactory,
                                              &finiteElementName]
                                               ( localIndex const targetRegionIndex, ElementSubRegionBase & elementSubRegion )
  {
    localIndex const numElems = elementSubRegion.size();
 
    // Get the constitutive model...and allocate a null constitutive model if required.
    constitutive::ConstitutiveBase * constitutiveRelation = nullptr;
    constitutive::NullModel * nullConstitutiveModel = nullptr;
    if( elementSubRegion.template hasWrapper< string >( constitutiveStringName ) )
    {
      string const & constitutiveName = elementSubRegion.template getReference< string >( constitutiveStringName );
      constitutiveRelation = &elementSubRegion.template getConstitutiveModel< >( constitutiveName );
    }
    else
    {
      nullConstitutiveModel = &elementSubRegion.template registerGroup< constitutive::NullModel >( "nullModelGroup" );
      constitutiveRelation = nullConstitutiveModel;
    }
 
    FiniteElementBase & subRegionFE = elementSubRegion.template getReference< FiniteElementBase >( finiteElementName );
 
    auto kernelLaunch = [&]( auto typeCombination )
    {
      using SUBREGION_TYPE = camp::at_v< decltype( typeCombination ), 0 >;
      using CONSTITUTIVE_TYPE = camp::at_v< decltype( typeCombination ), 1 >;
      using FE_TYPE = camp::at_v< decltype( typeCombination ), 2 >;
 
      SUBREGION_TYPE & castedSubRegion = dynamicCast< SUBREGION_TYPE & >( elementSubRegion );
      CONSTITUTIVE_TYPE & castedConstitutiveRelation = dynamicCast< CONSTITUTIVE_TYPE & >( *constitutiveRelation );
      typename FE_TYPE::ImplType & castedFiniteElement = *(dynamicCast< FE_TYPE & >( subRegionFE ).getImpl());
 
      auto kernel = kernelFactory.createKernel( nodeManager,
                                                edgeManager,
                                                faceManager,
                                                targetRegionIndex,
                                                castedSubRegion,
                                                castedFiniteElement,
                                                castedConstitutiveRelation );
 
      using KERNEL_TYPE = decltype( kernel );
 
      // Call the kernelLaunch function, and store the maximum contribution to the residual.
      maxResidualContribution =
        std::max( maxResidualContribution,
                  KERNEL_TYPE::template kernelLaunch< POLICY, KERNEL_TYPE >( numElems, kernel ) );
    };
 
    types::dispatch( DISPATCH_TYPE_LIST{}, kernelLaunch, elementSubRegion, *constitutiveRelation, subRegionFE );
 
    // Remove the null constitutive model (not required, but cleaner)
    if( nullConstitutiveModel )
    {
      elementSubRegion.deregisterGroup( "nullModelGroup" );
    }
 
  } );
 
  return maxResidualContribution;
}
//END_regionBasedKernelApplication
 
} // namespace finiteElement
} // namespace geos
 
 
 
#endif /* GEOS_FINITEELEMENT_KERNELBASE_HPP_ */