GEOS
PoromechanicsBase.hpp
Go to the documentation of this file.
1 /*
2  * ------------------------------------------------------------------------------------------------------------
3  * SPDX-License-Identifier: LGPL-2.1-only
4  *
5  * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
6  * Copyright (c) 2018-2024 TotalEnergies
7  * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
8  * Copyright (c) 2023-2024 Chevron
9  * Copyright (c) 2019- GEOS/GEOSX Contributors
10  * All rights reserved
11  *
12  * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
13  * ------------------------------------------------------------------------------------------------------------
14  */
15 
20 #ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_POROMECHANICSBASE_HPP_
21 #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_POROMECHANICSBASE_HPP_
22 
23 #include "finiteElement/BilinearFormUtilities.hpp"
24 #include "finiteElement/LinearFormUtilities.hpp"
28 
29 namespace geos
30 {
31 
32 namespace poromechanicsKernels
33 {
34 
50 template< typename SUBREGION_TYPE,
51  typename CONSTITUTIVE_TYPE,
52  typename FE_TYPE >
54  public finiteElement::ImplicitKernelBase< SUBREGION_TYPE,
55  CONSTITUTIVE_TYPE,
56  FE_TYPE,
57  3,
58  3 >
59 {
60 public:
61 
63  using Base = finiteElement::ImplicitKernelBase< SUBREGION_TYPE,
64  CONSTITUTIVE_TYPE,
65  FE_TYPE,
66  3,
67  3 >;
68 
75  using Base::m_dofNumber;
79  using Base::m_meshData;
80 
86  PoromechanicsBase( NodeManager const & nodeManager,
87  EdgeManager const & edgeManager,
88  FaceManager const & faceManager,
89  localIndex const targetRegionIndex,
90  SUBREGION_TYPE const & elementSubRegion,
91  FE_TYPE const & finiteElementSpace,
92  CONSTITUTIVE_TYPE & inputConstitutiveType,
93  arrayView1d< globalIndex const > const inputDispDofNumber,
94  globalIndex const rankOffset,
96  arrayView1d< real64 > const inputRhs,
97  real64 const inputDt,
98  real64 const (&gravityVector)[3],
99  string const inputFlowDofKey,
100  string const fluidModelKey ):
101  Base( nodeManager,
102  edgeManager,
103  faceManager,
104  targetRegionIndex,
105  elementSubRegion,
106  finiteElementSpace,
107  inputConstitutiveType,
108  inputDispDofNumber,
109  rankOffset,
110  inputMatrix,
111  inputRhs,
112  inputDt ),
113  m_X( nodeManager.referencePosition() ),
114  m_disp( nodeManager.getField< fields::solidMechanics::totalDisplacement >() ),
115  m_uhat( nodeManager.getField< fields::solidMechanics::incrementalDisplacement >() ),
116  m_gravityVector{ gravityVector[0], gravityVector[1], gravityVector[2] },
117  m_gravityAcceleration( LvArray::tensorOps::l2Norm< 3 >( gravityVector ) ),
118  m_solidDensity( inputConstitutiveType.getDensity() ),
119  m_flowDofNumber( elementSubRegion.template getReference< array1d< globalIndex > >( inputFlowDofKey ) ),
120  m_pressure_n( elementSubRegion.template getField< fields::flow::pressure_n >() ),
121  m_pressure( elementSubRegion.template getField< fields::flow::pressure >() )
122  {
123  GEOS_UNUSED_VAR( fluidModelKey );
124  }
125 
126  //*****************************************************************************
135  {
136 public:
137 
138  static constexpr int numDispDofPerElem = Base::StackVariables::maxNumRows;
139 
143  Base::StackVariables(),
144  xLocal(),
145  u_local(),
146  uhat_local(),
147  localResidualMomentum( Base::StackVariables::localResidual ),
150  localPressureDofIndex{ 0 }
151  {}
152 
153 #if !defined(CALC_FEM_SHAPE_IN_KERNEL)
155  int xLocal;
156 #else
158  real64 xLocal[numNodesPerElem][3];
159 #endif
160 
161  // Storage for displacements
162 
167 
168 
169  // Storage for helper variables used in the quadrature point kernel
170 
172  real64 strainIncrement[6]{};
173 
175  real64 totalStress[6]{};
177  typename CONSTITUTIVE_TYPE::KernelWrapper::DiscretizationOps stiffness;
179  real64 dTotalStress_dPressure[6]{};
181  real64 dTotalStress_dTemperature[6]{};
182 
184  real64 bodyForce[3]{};
186  real64 dBodyForce_dVolStrainIncrement[3]{};
188  real64 dBodyForce_dPressure[3]{};
189 
191  real64 temperature{}; // for stress computation
193  real64 deltaTemperatureFromLastStep{}; // for porosity update
194 
195  // Storage for residual and degrees of freedom
196 
198  real64 ( &localResidualMomentum )[numDispDofPerElem];
200  real64 ( &dLocalResidualMomentum_dDisplacement )[numDispDofPerElem][numDispDofPerElem];
202  real64 dLocalResidualMomentum_dPressure[numDispDofPerElem][1]{};
203 
205  globalIndex localPressureDofIndex{};
206 
207  };
208  //*****************************************************************************
209 
220  void setup( localIndex const k,
221  StackVariables & stack ) const
222  {
223  integer constexpr numDims = 3;
224 
225  m_finiteElementSpace.template setup< FE_TYPE >( k, m_meshData, stack.feStack );
226  localIndex const numSupportPoints =
227  m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
228 
229  for( localIndex a=0; a<numSupportPoints; ++a )
230  {
231  localIndex const localNodeIndex = m_elemsToNodes( k, a );
232 
233  for( integer i = 0; i < numDims; ++i )
234  {
235 #if defined(CALC_FEM_SHAPE_IN_KERNEL)
236  stack.xLocal[a][i] = m_X[localNodeIndex][i];
237 #endif
238  stack.u_local[a][i] = m_disp[localNodeIndex][i];
239  stack.uhat_local[a][i] = m_uhat[localNodeIndex][i];
240  stack.localRowDofIndex[a*numDims+i] = m_dofNumber[localNodeIndex]+i;
241  stack.localColDofIndex[a*numDims+i] = m_dofNumber[localNodeIndex]+i;
242  }
243  }
244 
245  stack.localPressureDofIndex = m_flowDofNumber[k];
246 
247  // Add stabilization to block diagonal parts of the local dResidualMomentum_dDisplacement (this
248  // is a no-operation with FEM classes)
249  real64 const stabilizationScaling = computeStabilizationScaling( k );
250  m_finiteElementSpace.template addGradGradStabilizationMatrix
251  < FE_TYPE, numDofPerTrialSupportPoint, false >( stack.feStack,
253  -stabilizationScaling );
254  m_finiteElementSpace.template
255  addEvaluatedGradGradStabilizationVector< FE_TYPE,
256  numDofPerTrialSupportPoint >
257  ( stack.feStack,
258  stack.uhat_local,
259  reinterpret_cast< real64 (&)[numNodesPerElem][numDofPerTestSupportPoint] >(stack.localResidualMomentum),
260  -stabilizationScaling );
261  }
262 
272  {
273  // TODO: generalize this to other constitutive models (currently we assume linear elasticity).
274  return 2.0 * m_constitutiveUpdate.getShearModulus( k );
275  }
276 
277 protected:
278 
281 
284 
287 
289  real64 const m_gravityVector[3]{};
292 
295 
298 
303 
304 };
305 
306 
307 } // namespace poromechanicsKernels
308 
309 } // namespace geos
310 
311 #endif // GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_POROMECHANICSBASE_HPP_
#define GEOS_HOST_DEVICE
Marks a host-device function.
Definition: GeosxMacros.hpp:49
#define GEOS_UNUSED_VAR(...)
Mark an unused variable and silence compiler warnings.
Definition: GeosxMacros.hpp:84
#define GEOS_FORCE_INLINE
Marks a function or lambda for inlining.
Definition: GeosxMacros.hpp:51
This class provides an interface to ObjectManagerBase in order to manage edge data.
Definition: EdgeManager.hpp:43
The FaceManager class provides an interface to ObjectManagerBase in order to manage face data.
Definition: FaceManager.hpp:44
The NodeManager class provides an interface to ObjectManagerBase in order to manage node data.
Definition: NodeManager.hpp:46
Define the base interface for implicit finite element kernels.
traits::ViewTypeConst< typename SUBREGION_TYPE::NodeMapType::base_type > const m_elemsToNodes
The element to nodes map.
Definition: KernelBase.hpp:258
FE_TYPE::template MeshData< SUBREGION_TYPE > m_meshData
Data structure containing mesh data used to setup the finite element.
arrayView1d< globalIndex const > const m_dofNumber
The global degree of freedom number.
CONSTITUTIVE_TYPE::KernelWrapper const m_constitutiveUpdate
Definition: KernelBase.hpp:265
Defines the kernel structure for solving quasi-static poromechanics.
arrayView1d< real64 const > const m_pressure
The rank-global fluid pressure.
arrayView2d< real64 const > m_solidDensity
The rank global density.
GEOS_HOST_DEVICE GEOS_FORCE_INLINE void setup(localIndex const k, StackVariables &stack) const
Copy global values from primary field to a local stack array.
arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const m_X
The array containing the nodal position array.
arrayView2d< real64 const, nodes::INCR_DISPLACEMENT_USD > const m_uhat
The rank-global incremental displacement array.
real64 const m_gravityVector[3]
The gravity vector.
arrayView2d< real64 const, nodes::TOTAL_DISPLACEMENT_USD > const m_disp
The rank-global displacement array.
GEOS_HOST_DEVICE GEOS_FORCE_INLINE real64 computeStabilizationScaling(localIndex const k) const
Get a parameter representative of the stiffness, used as physical scaling for the stabilization matri...
arrayView1d< globalIndex const > const m_flowDofNumber
The global degree of freedom number.
real64 const m_gravityAcceleration
The L2-norm of the gravity vector.
arrayView1d< real64 const > const m_pressure_n
The rank-global fluid pressure at the previous converged time step.
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)
Constructor.
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
Definition: DataTypes.hpp:188
LvArray::CRSMatrixView< T, COL_INDEX, localIndex const, LvArray::ChaiBuffer > CRSMatrixView
Alias for CRS Matrix View.
Definition: DataTypes.hpp:318
GEOS_GLOBALINDEX_TYPE globalIndex
Global index type (for indexing objects across MPI partitions).
Definition: DataTypes.hpp:87
double real64
64-bit floating point type.
Definition: DataTypes.hpp:98
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Definition: DataTypes.hpp:84
std::int32_t integer
Signed integer type.
Definition: DataTypes.hpp:81
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
Definition: DataTypes.hpp:204
Array< T, 1 > array1d
Alias for 1D array.
Definition: DataTypes.hpp:184
Kernel variables allocated on the stack.
Definition: KernelBase.hpp:138
real64(& dLocalResidualMomentum_dDisplacement)[numDispDofPerElem][numDispDofPerElem]
Derivative of linear momentum balance residual wrt displacement.
globalIndex localPressureDofIndex
C-array storage for the element local row degrees of freedom.
CONSTITUTIVE_TYPE::KernelWrapper::DiscretizationOps stiffness
Derivative of total stress wrt displacement.
real64 uhat_local[numNodesPerElem][numDofPerTrialSupportPoint]
Stack storage for the element local nodal incremental displacement.
real64 u_local[numNodesPerElem][numDofPerTrialSupportPoint]
Stack storage for the element local nodal displacement.
real64 dLocalResidualMomentum_dPressure[numDispDofPerElem][1]
Derivative of linear momentum balance residual wrt pressure.
real64(& localResidualMomentum)[numDispDofPerElem]
Linear momentum balance residual.