ImplicitKernelBase.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.
 * ------------------------------------------------------------------------------------------------------------
 */
 
#include "KernelBase.hpp"
 
#ifndef GEOS_FINITEELEMENT_IMPLICITKERNELBASE_HPP_
#define GEOS_FINITEELEMENT_IMPLICITKERNELBASE_HPP_
 
 
 
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 ImplicitKernelBase : public KernelBase< SUBREGION_TYPE,
                                              CONSTITUTIVE_TYPE,
                                              FE_TYPE,
                                              NUM_DOF_PER_TEST_SP,
                                              NUM_DOF_PER_TRIAL_SP >
{
public:
  using Base = KernelBase< SUBREGION_TYPE,
                           CONSTITUTIVE_TYPE,
                           FE_TYPE,
                           NUM_DOF_PER_TEST_SP,
                           NUM_DOF_PER_TRIAL_SP >;
 
  using Base::maxNumTestSupportPointsPerElem;
  using Base::maxNumTrialSupportPointsPerElem;
  using Base::numDofPerTestSupportPoint;
  using Base::numDofPerTrialSupportPoint;
  using Base::m_elemsToNodes;
  using Base::m_finiteElementSpace;
 
 
  ImplicitKernelBase( NodeManager const & nodeManager,
                      EdgeManager const & edgeManager,
                      FaceManager const & faceManager,
                      localIndex const targetRegionIndex,
                      SUBREGION_TYPE const & elementSubRegion,
                      FE_TYPE const & finiteElementSpace,
                      CONSTITUTIVE_TYPE & inputConstitutiveType,
                      arrayView1d< globalIndex const > const & inputDofNumber,
                      globalIndex const rankOffset,
                      CRSMatrixView< real64, globalIndex const > const & inputMatrix,
                      arrayView1d< real64 > const & inputRhs,
                      real64 const inputDt ):
    Base( elementSubRegion,
          finiteElementSpace,
          inputConstitutiveType ),
    m_dofNumber( inputDofNumber ),
    m_dofRankOffset( rankOffset ),
    m_matrix( inputMatrix ),
    m_rhs( inputRhs ),
    m_dt( inputDt )
  {
    FiniteElementBase::initialize< FE_TYPE >( nodeManager,
                                              edgeManager,
                                              faceManager,
                                              elementSubRegion,
                                              m_meshData );
    GEOS_UNUSED_VAR( targetRegionIndex );
  }
 
 
  //***************************************************************************
  struct StackVariables : public Base::StackVariables
  {
    static constexpr int maxNumRows = maxNumTestSupportPointsPerElem *numDofPerTestSupportPoint;
 
    static constexpr int maxNumCols = maxNumTrialSupportPointsPerElem *numDofPerTrialSupportPoint;
 
    GEOS_HOST_DEVICE
    StackVariables():
      localRowDofIndex{ 0 },
      localColDofIndex{ 0 },
      localResidual{ 0.0 },
      localJacobian{ {0.0} }
    {
      for( int ii = 0; ii < maxNumRows; ++ii )
      {
        for( int jj = 0; jj < maxNumCols; ++jj )
        {
          localJacobian[ii][jj] = 0.0;
        }
      }
    }
 
    localIndex numRows;
 
    localIndex numCols;
 
    globalIndex localRowDofIndex[maxNumRows];
 
    globalIndex localColDofIndex[maxNumCols];
 
    real64 localResidual[maxNumRows];
 
    real64 localJacobian[maxNumRows][maxNumCols];
 
    typename FE_TYPE::StackVariables feStack;
  };
  //***************************************************************************
 
  GEOS_HOST_DEVICE
  inline
  void setup( localIndex const k,
              StackVariables & stack ) const
  {
    m_finiteElementSpace.template setup< FE_TYPE >( k, m_meshData, stack.feStack );
    localIndex numTestSupportPoints = m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
    localIndex numTrialSupportPoints = numTestSupportPoints;
    stack.numRows = numTestSupportPoints * numDofPerTestSupportPoint;
    stack.numCols = numTrialSupportPoints * numDofPerTrialSupportPoint;
    for( localIndex a=0; a<numTestSupportPoints; ++a )
    {
      localIndex const localNodeIndex = m_elemsToNodes[k][a];
      for( int i=0; i<numDofPerTestSupportPoint; ++i )
      {
        stack.localRowDofIndex[a*numDofPerTestSupportPoint+i] = m_dofNumber[localNodeIndex]+i;
      }
    }
 
    for( localIndex a=0; a<numTrialSupportPoints; ++a )
    {
      localIndex const localNodeIndex = m_elemsToNodes[k][a];
      for( int i=0; i<numDofPerTrialSupportPoint; ++i )
      {
        stack.localColDofIndex[a*numDofPerTrialSupportPoint+i] = m_dofNumber[localNodeIndex]+i;
      }
    }
  }
 
 
protected:
  arrayView1d< globalIndex const > const m_dofNumber;
 
  globalIndex const m_dofRankOffset;
 
  CRSMatrixView< real64, globalIndex const > const m_matrix;
 
  arrayView1d< real64 > const m_rhs;
 
  typename FE_TYPE::template MeshData< SUBREGION_TYPE > m_meshData;
 
  real64 const m_dt; 
 
};
 
}
}
 
 
 
#endif /* GEOS_FINITEELEMENT_IMPLICITKERNELBASE_HPP_ */