PoromechanicsBase.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_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_POROMECHANICSBASE_HPP_
#define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_POROMECHANICSBASE_HPP_
 
#include "finiteElement/BilinearFormUtilities.hpp"
#include "finiteElement/LinearFormUtilities.hpp"
#include "finiteElement/kernelInterface/ImplicitKernelBase.hpp"
#include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp"
#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp"
 
namespace geos
{
 
namespace poromechanicsKernels
{
 
template< typename SUBREGION_TYPE,
          typename CONSTITUTIVE_TYPE,
          typename FE_TYPE >
class PoromechanicsBase :
  public finiteElement::ImplicitKernelBase< SUBREGION_TYPE,
                                            CONSTITUTIVE_TYPE,
                                            FE_TYPE,
                                            3,
                                            3 >
{
public:
 
  using Base = finiteElement::ImplicitKernelBase< SUBREGION_TYPE,
                                                  CONSTITUTIVE_TYPE,
                                                  FE_TYPE,
                                                  3,
                                                  3 >;
 
  static constexpr int numNodesPerElem = Base::maxNumTestSupportPointsPerElem;
  using Base::numDofPerTestSupportPoint;
  using Base::numDofPerTrialSupportPoint;
  using Base::m_dofNumber;
  using Base::m_elemsToNodes;
  using Base::m_constitutiveUpdate;
  using Base::m_finiteElementSpace;
  using Base::m_meshData;
 
  PoromechanicsBase( 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 inputDispDofNumber,
                     globalIndex const rankOffset,
                     CRSMatrixView< real64, globalIndex const > const inputMatrix,
                     arrayView1d< real64 > const inputRhs,
                     real64 const inputDt,
                     real64 const (&gravityVector)[3],
                     string const inputFlowDofKey,
                     string const fluidModelKey ):
    Base( nodeManager,
          edgeManager,
          faceManager,
          targetRegionIndex,
          elementSubRegion,
          finiteElementSpace,
          inputConstitutiveType,
          inputDispDofNumber,
          rankOffset,
          inputMatrix,
          inputRhs,
          inputDt ),
    m_X( nodeManager.referencePosition() ),
    m_disp( nodeManager.getField< fields::solidMechanics::totalDisplacement >() ),
    m_uhat( nodeManager.getField< fields::solidMechanics::incrementalDisplacement >() ),
    m_gravityVector{ gravityVector[0], gravityVector[1], gravityVector[2] },
    m_gravityAcceleration( LvArray::tensorOps::l2Norm< 3 >( gravityVector ) ),
    m_solidDensity( inputConstitutiveType.getDensity() ),
    m_flowDofNumber( elementSubRegion.template getReference< array1d< globalIndex > >( inputFlowDofKey ) ),
    m_pressure_n( elementSubRegion.template getField< fields::flow::pressure_n >() ),
    m_pressure( elementSubRegion.template getField< fields::flow::pressure >() )
  {
    GEOS_UNUSED_VAR( fluidModelKey );
  }
 
  //*****************************************************************************
  struct StackVariables : public Base::StackVariables
  {
public:
 
    static constexpr int numDispDofPerElem =  Base::StackVariables::maxNumRows;
 
    GEOS_HOST_DEVICE
    StackVariables():
      Base::StackVariables(),
            xLocal(),
            u_local(),
            uhat_local(),
            localResidualMomentum( Base::StackVariables::localResidual ),
      dLocalResidualMomentum_dDisplacement( Base::StackVariables::localJacobian ),
      dLocalResidualMomentum_dPressure{ {0.0} },
      localPressureDofIndex{ 0 }
    {}
 
#if !defined(CALC_FEM_SHAPE_IN_KERNEL)
    int xLocal;
#else
    real64 xLocal[numNodesPerElem][3];
#endif
 
    // Storage for displacements
 
    real64 u_local[numNodesPerElem][numDofPerTrialSupportPoint]{};
    real64 uhat_local[numNodesPerElem][numDofPerTrialSupportPoint]{};
 
 
    // Storage for helper variables used in the quadrature point kernel
 
    real64 strainIncrement[6]{};
 
    real64 totalStress[6]{};
    typename CONSTITUTIVE_TYPE::KernelWrapper::DiscretizationOps stiffness;
    real64 dTotalStress_dPressure[6]{};
    real64 dTotalStress_dTemperature[6]{};
 
    real64 bodyForce[3]{};
    real64 dBodyForce_dVolStrainIncrement[3]{};
    real64 dBodyForce_dPressure[3]{};
 
    real64 temperature{}; // for stress computation
    real64 deltaTemperatureFromLastStep{}; // for porosity update
    real64 deltaTemperature{}; // for stress computation
 
 
    // Storage for residual and degrees of freedom
 
    real64 ( &localResidualMomentum )[numDispDofPerElem];
    real64 ( &dLocalResidualMomentum_dDisplacement )[numDispDofPerElem][numDispDofPerElem];
    real64 dLocalResidualMomentum_dPressure[numDispDofPerElem][1]{};
 
    globalIndex localPressureDofIndex{};
 
  };
  //*****************************************************************************
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  void setup( localIndex const k,
              StackVariables & stack ) const
  {
    integer constexpr numDims = 3;
 
    m_finiteElementSpace.template setup< FE_TYPE >( k, m_meshData, stack.feStack );
    localIndex const numSupportPoints =
      m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
 
    for( localIndex a=0; a<numSupportPoints; ++a )
    {
      localIndex const localNodeIndex = m_elemsToNodes( k, a );
 
      for( integer i = 0; i < numDims; ++i )
      {
#if defined(CALC_FEM_SHAPE_IN_KERNEL)
        stack.xLocal[a][i] = m_X[localNodeIndex][i];
#endif
        stack.u_local[a][i] = m_disp[localNodeIndex][i];
        stack.uhat_local[a][i] = m_uhat[localNodeIndex][i];
        stack.localRowDofIndex[a*numDims+i] = m_dofNumber[localNodeIndex]+i;
        stack.localColDofIndex[a*numDims+i] = m_dofNumber[localNodeIndex]+i;
      }
    }
 
    stack.localPressureDofIndex = m_flowDofNumber[k];
 
    // Add stabilization to block diagonal parts of the local dResidualMomentum_dDisplacement (this
    // is a no-operation with FEM classes)
    real64 const stabilizationScaling = computeStabilizationScaling( k );
    m_finiteElementSpace.template addGradGradStabilizationMatrix
    < FE_TYPE, numDofPerTrialSupportPoint, false >( stack.feStack,
                                                    stack.dLocalResidualMomentum_dDisplacement,
                                                    -stabilizationScaling );
    m_finiteElementSpace.template
    addEvaluatedGradGradStabilizationVector< FE_TYPE,
                                             numDofPerTrialSupportPoint >
      ( stack.feStack,
      stack.uhat_local,
      reinterpret_cast< real64 (&)[numNodesPerElem][numDofPerTestSupportPoint] >(stack.localResidualMomentum),
      -stabilizationScaling );
  }
 
  GEOS_HOST_DEVICE
  GEOS_FORCE_INLINE
  real64 computeStabilizationScaling( localIndex const k ) const
  {
    // TODO: generalize this to other constitutive models (currently we assume linear elasticity).
    return 2.0 * m_constitutiveUpdate.getShearModulus( k );
  }
 
protected:
 
  arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const m_X;
 
  arrayView2d< real64 const, nodes::TOTAL_DISPLACEMENT_USD > const m_disp;
 
  arrayView2d< real64 const, nodes::INCR_DISPLACEMENT_USD > const m_uhat;
 
  real64 const m_gravityVector[3]{};
  real64 const m_gravityAcceleration;
 
  arrayView2d< real64 const > m_solidDensity;
 
  arrayView1d< globalIndex const > const m_flowDofNumber;
 
  arrayView1d< real64 const > const m_pressure_n;
  arrayView1d< real64 const > const m_pressure;
 
};
 
 
} // namespace poromechanicsKernels
 
} // namespace geos
 
#endif // GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_POROMECHANICSBASE_HPP_