KernelBase.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 Total, S.A
 * Copyright (c) 2019-     GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */


#ifndef GEOSX_FINITEELEMENT_KERNELBASE_HPP_
#define GEOSX_FINITEELEMENT_KERNELBASE_HPP_

#include "common/DataTypes.hpp"
#include "common/TimingMacros.hpp"
#include "constitutive/ConstitutivePassThru.hpp"
#include "finiteElement/FiniteElementDispatch.hpp"
#include "mesh/ElementRegionManager.hpp"
#include "rajaInterface/GEOS_RAJA_Interface.hpp"



#if defined(__APPLE__)
#define CONSTRUCTOR_PARAM_OPTION 2
#else
#define CONSTRUCTOR_PARAM_OPTION 1
#endif

#if CONSTRUCTOR_PARAM_OPTION==1
namespace std
{
namespace detail
{
template< class T, class Tuple, std::size_t... I >
constexpr T make_from_tuple_impl( Tuple && t, std::index_sequence< I... > )
{
  return T( std::get< I >( std::forward< Tuple >( t ))... );
}
} // namespace detail

template< class T, class Tuple >
constexpr T make_from_tuple( Tuple && t )
{
  return detail::make_from_tuple_impl< T >( std::forward< Tuple >( t ),
                                            std::make_index_sequence< std::tuple_size< std::remove_reference_t< Tuple > >::value >{} );
}

}
#elif CONSTRUCTOR_PARAM_OPTION==2
  #include "camp/camp.hpp"
#endif

namespace geosx
{

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 numTestSupportPointsPerElem  = FE_TYPE::numNodes;

  static constexpr int numTrialSupportPointsPerElem = FE_TYPE::numNodes;

  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 * const inputConstitutiveType ):
    m_elemsToNodes( elementSubRegion.nodeList().toViewConst() ),
    m_elemGhostRank( elementSubRegion.ghostRank() ),
    m_constitutiveUpdate( inputConstitutiveType->createKernelUpdates() ),
    m_finiteElementSpace( finiteElementSpace )
  {}

  struct StackVariables
  {};

  GEOSX_HOST_DEVICE
  GEOSX_FORCE_INLINE
  void setup( localIndex const k,
              StackVariables & stack ) const
  {
    GEOSX_UNUSED_VAR( k );
    GEOSX_UNUSED_VAR( stack );
  }

  GEOSX_HOST_DEVICE
  GEOSX_FORCE_INLINE
  void quadraturePointKernel( localIndex const k,
                              localIndex const q,
                              StackVariables & stack ) const
  {
    GEOSX_UNUSED_VAR( k );
    GEOSX_UNUSED_VAR( q );
    GEOSX_UNUSED_VAR( stack );
  }

  GEOSX_HOST_DEVICE
  GEOSX_FORCE_INLINE
  real64 complete( localIndex const k,
                   StackVariables & stack ) const
  {
    GEOSX_UNUSED_VAR( k );
    GEOSX_UNUSED_VAR( stack );
    return 0;
  }


  //START_kernelLauncher
  template< typename POLICY,
            typename KERNEL_TYPE >
  static
  real64
  kernelLaunch( localIndex const numElems,
                KERNEL_TYPE const & kernelComponent )
  {
    GEOSX_MARK_FUNCTION;

    // Define a RAJA reduction variable to get the maximum residual contribution.
    RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxResidual( 0 );

    forAll< POLICY >( numElems,
                      [=] GEOSX_HOST_DEVICE ( localIndex const k )
    {
      typename KERNEL_TYPE::StackVariables stack;

      kernelComponent.setup( k, stack );
      for( integer q=0; q<numQuadraturePointsPerElem; ++q )
      {
        kernelComponent.quadraturePointKernel( k, q, stack );
      }
      maxResidual.max( kernelComponent.complete( k, stack ) );
    } );
    return maxResidual.get();
  }
  //END_kernelLauncher

protected:
  traits::ViewTypeConst< typename SUBREGION_TYPE::NodeMapType::base_type > const m_elemsToNodes;

  arrayView1d< integer const > const m_elemGhostRank;

  typename CONSTITUTIVE_TYPE::KernelWrapper const m_constitutiveUpdate;

  FE_TYPE const & m_finiteElementSpace;
};


//*****************************************************************************
//*****************************************************************************
//*****************************************************************************

//START_regionBasedKernelApplication
template< typename POLICY,
          typename CONSTITUTIVE_BASE,
          typename REGION_TYPE,
          template< typename SUBREGION_TYPE,
                    typename CONSTITUTIVE_TYPE,
                    typename FE_TYPE > class KERNEL_TEMPLATE,
          typename ... KERNEL_CONSTRUCTOR_PARAMS >
static
real64 regionBasedKernelApplication( MeshLevel & mesh,
                                     arrayView1d< string const > const & targetRegions,
                                     string const & finiteElementName,
                                     arrayView1d< string const > const & constitutiveNames,
                                     KERNEL_CONSTRUCTOR_PARAMS && ... kernelConstructorParams )
{
  GEOSX_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());


  // Create a tuple that contains the kernelConstructorParams, as the lambda does not properly catch the parameter pack
  // until c++20
#if CONSTRUCTOR_PARAM_OPTION==1
  std::tuple< KERNEL_CONSTRUCTOR_PARAMS &... > kernelConstructorParamsTuple = std::forward_as_tuple( kernelConstructorParams ... );
#elif CONSTRUCTOR_PARAM_OPTION==2
  camp::tuple< KERNEL_CONSTRUCTOR_PARAMS &... > kernelConstructorParamsTuple = camp::forward_as_tuple( kernelConstructorParams ... );
#endif


  // Loop over all sub-regions in regiongs of type REGION_TYPE, that are listed in the targetRegions array.
  elementRegionManager.forElementSubRegions< REGION_TYPE >( targetRegions,
                                                            [&constitutiveNames,
                                                             &maxResidualContribution,
                                                             &nodeManager,
                                                             &edgeManager,
                                                             &faceManager,
                                                             &kernelConstructorParamsTuple,
                                                             &finiteElementName]
                                                              ( localIndex const targetRegionIndex, auto & elementSubRegion )
  {
    localIndex const numElems = elementSubRegion.size();

    // Create an alias for the type of subregion we are in, which is now known at compile time.
    typedef TYPEOFREF( elementSubRegion ) SUBREGIONTYPE;

    // Get the constitutive model...and allocate a null constitutive model if required.
    constitutive::ConstitutiveBase * constitutiveRelation = nullptr;
    constitutive::NullModel * nullConstitutiveModel = nullptr;
    if( targetRegionIndex <= constitutiveNames.size()-1 )
    {
      constitutiveRelation = elementSubRegion.template getConstitutiveModel( constitutiveNames[targetRegionIndex] );
    }
    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,
                                                                       &kernelConstructorParamsTuple,
                                                                       &elementSubRegion,
                                                                       &finiteElementName,
                                                                       numElems]
                                                                        ( auto * const castedConstitutiveRelation )
    {
      // Create an alias for the type of constitutive model.
      using CONSTITUTIVE_TYPE = TYPEOFPTR( castedConstitutiveRelation );


      string const elementTypeString = elementSubRegion.GetElementTypeString();

      FiniteElementBase &
      subRegionFE = elementSubRegion.template getReference< FiniteElementBase >( finiteElementName );

      finiteElement::dispatch3D( subRegionFE,
                                 [&maxResidualContribution,
                                  &nodeManager,
                                  &edgeManager,
                                  &faceManager,
                                  &kernelConstructorParamsTuple,
                                  &elementSubRegion,
                                  &numElems,
                                  &castedConstitutiveRelation] ( auto const finiteElement )
      {
        using FE_TYPE = TYPEOFREF( finiteElement );

        // Define an alias for the kernel type for easy use.
        using KERNEL_TYPE = KERNEL_TEMPLATE< SUBREGIONTYPE,
                                             CONSTITUTIVE_TYPE,
                                             FE_TYPE >;

        // 1) Combine the tuple containing the physics kernel specific constructor parameters with
        // the parameters common to all phsyics kernels that use this interface,
        // 2) Instantiate the kernel.
        // note: have two options, using std::tuple and camp::tuple. Due to a bug in the OSX
        // implementation of std::tuple_cat, we must use camp on OSX. In the future, we should
        // only use one option...most likely camp since we can easily fix bugs.
#if CONSTRUCTOR_PARAM_OPTION==1
        auto temp = std::forward_as_tuple( nodeManager,
                                           edgeManager,
                                           faceManager,
                                           elementSubRegion,
                                           finiteElement,
                                           castedConstitutiveRelation );

        auto fullKernelComponentConstructorArgs = std::tuple_cat( temp,
                                                                  kernelConstructorParamsTuple );

        KERNEL_TYPE kernelComponent = std::make_from_tuple< KERNEL_TYPE >( fullKernelComponentConstructorArgs );

#elif CONSTRUCTOR_PARAM_OPTION==2
        auto temp = camp::forward_as_tuple( nodeManager,
                                            edgeManager,
                                            faceManager,
                                            elementSubRegion,
                                            finiteElement,
                                            castedConstitutiveRelation );
        auto fullKernelComponentConstructorArgs = camp::tuple_cat_pair( temp,
                                                                        kernelConstructorParamsTuple );
        KERNEL_TYPE kernelComponent = camp::make_from_tuple< KERNEL_TYPE >( fullKernelComponentConstructorArgs );

#endif

        // Call the kernelLaunch function, and store the maximum contribution to the residual.
        maxResidualContribution =
          std::max( maxResidualContribution,
                    KERNEL_TYPE::template kernelLaunch< POLICY,
                                                        KERNEL_TYPE >( numElems,
                                                                       kernelComponent ) );
      } );
    } );

    // Remove the null constitutive model (not required, but cleaner)
    if( nullConstitutiveModel )
    {
      elementSubRegion.deregisterGroup( "nullModelGroup" );
    }

  } );

  return maxResidualContribution;
}
//END_regionBasedKernelApplication

} // namespace finiteElement
} // namespace geosx



#endif /* GEOSX_FINITEELEMENT_KERNELBASE_HPP_ */