GEOS
SinglePhasePoromechanics_impl.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_SINGLEPHASEPOROMECHANICS_IMPL_HPP_
21 #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_SINGLEPHASEPOROMECHANICS_IMPL_HPP_
22 
23 #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp"
24 #include "finiteElement/BilinearFormUtilities.hpp"
25 #include "finiteElement/LinearFormUtilities.hpp"
27 
28 namespace geos
29 {
30 
31 namespace poromechanicsKernels
32 {
33 
34 template< typename SUBREGION_TYPE,
35  typename CONSTITUTIVE_TYPE,
36  typename FE_TYPE >
38 SinglePhasePoromechanics( NodeManager const & nodeManager,
39  EdgeManager const & edgeManager,
40  FaceManager const & faceManager,
41  localIndex const targetRegionIndex,
42  SUBREGION_TYPE const & elementSubRegion,
43  FE_TYPE const & finiteElementSpace,
44  CONSTITUTIVE_TYPE & inputConstitutiveType,
45  arrayView1d< globalIndex const > const inputDispDofNumber,
46  globalIndex const rankOffset,
48  arrayView1d< real64 > const inputRhs,
49  real64 const inputDt,
50  real64 const (&gravityVector)[3],
51  string const inputFlowDofKey,
52  integer const performStressInitialization,
53  string const fluidModelKey ):
54  Base( nodeManager,
55  edgeManager,
56  faceManager,
57  targetRegionIndex,
58  elementSubRegion,
59  finiteElementSpace,
60  inputConstitutiveType,
61  inputDispDofNumber,
62  rankOffset,
63  inputMatrix,
64  inputRhs,
65  inputDt,
66  gravityVector,
67  inputFlowDofKey,
68  fluidModelKey ),
69  m_fluidDensity( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).density() ),
70  m_fluidDensity_n( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).density_n() ),
71  m_dFluidDensity( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).dDensity() ),
72  m_performStressInitialization( performStressInitialization )
73 {}
74 
75 template< typename SUBREGION_TYPE,
76  typename CONSTITUTIVE_TYPE,
77  typename FE_TYPE >
82  localIndex const q,
83  StackVariables & stack ) const
84 {
85  real64 porosity = 0.0;
86  real64 porosity_n = 0.0;
87  real64 dPorosity_dVolStrain = 0.0;
88  real64 dPorosity_dPressure = 0.0;
89  real64 dPorosity_dTemperature = 0.0;
90  real64 dSolidDensity_dPressure = 0.0;
91 
92  // Step 1: call the constitutive model to evaluate the total stress and compute porosity
93  m_constitutiveUpdate.smallStrainUpdatePoromechanics( k, q,
94  m_dt,
95  m_pressure[k],
96  m_pressure_n[k],
97  stack.temperature,
98  stack.deltaTemperatureFromLastStep,
99  stack.strainIncrement,
100  stack.totalStress,
101  stack.dTotalStress_dPressure,
102  stack.dTotalStress_dTemperature,
103  stack.stiffness,
104  m_performStressInitialization,
105  porosity,
106  porosity_n,
107  dPorosity_dVolStrain,
108  dPorosity_dPressure,
109  dPorosity_dTemperature,
110  dSolidDensity_dPressure );
111 
112  // Step 2: compute the body force
113  computeBodyForce( k, q,
114  porosity,
115  dPorosity_dVolStrain,
116  dPorosity_dPressure,
117  dPorosity_dTemperature,
118  dSolidDensity_dPressure,
119  stack );
120 
121  // Step 3: compute fluid mass increment
122  computeFluidIncrement( k, q,
123  porosity,
124  porosity_n,
125  dPorosity_dVolStrain,
126  dPorosity_dPressure,
127  dPorosity_dTemperature,
128  stack );
129 }
130 
131 template< typename SUBREGION_TYPE,
132  typename CONSTITUTIVE_TYPE,
133  typename FE_TYPE >
138  localIndex const q,
139  real64 const & porosity,
140  real64 const & dPorosity_dVolStrain,
141  real64 const & dPorosity_dPressure,
142  real64 const & dPorosity_dTemperature,
143  real64 const & dSolidDensity_dPressure,
144  StackVariables & stack ) const
145 {
146  GEOS_UNUSED_VAR( dPorosity_dTemperature );
147 
148  real64 const mixtureDensity = ( 1.0 - porosity ) * m_solidDensity( k, q ) + porosity * m_fluidDensity( k, q );
149  real64 const dMixtureDens_dVolStrainIncrement = dPorosity_dVolStrain * ( -m_solidDensity( k, q ) + m_fluidDensity( k, q ) );
150  real64 const dMixtureDens_dPressure = dPorosity_dPressure * ( -m_solidDensity( k, q ) + m_fluidDensity( k, q ) )
151  + ( 1.0 - porosity ) * dSolidDensity_dPressure
152  + porosity * m_dFluidDensity( k, q, DerivOffset::dP );
153 
154  LvArray::tensorOps::scaledCopy< 3 >( stack.bodyForce, m_gravityVector, mixtureDensity );
155  LvArray::tensorOps::scaledCopy< 3 >( stack.dBodyForce_dVolStrainIncrement, m_gravityVector, dMixtureDens_dVolStrainIncrement );
156  LvArray::tensorOps::scaledCopy< 3 >( stack.dBodyForce_dPressure, m_gravityVector, dMixtureDens_dPressure );
157 }
158 
159 template< typename SUBREGION_TYPE,
160  typename CONSTITUTIVE_TYPE,
161  typename FE_TYPE >
166  localIndex const q,
167  real64 const & porosity,
168  real64 const & porosity_n,
169  real64 const & dPorosity_dVolStrain,
170  real64 const & dPorosity_dPressure,
171  real64 const & dPorosity_dTemperature,
172  StackVariables & stack ) const
173 {
174  GEOS_UNUSED_VAR( dPorosity_dTemperature );
175 
176  stack.fluidMassIncrement = porosity * m_fluidDensity( k, q ) - porosity_n * m_fluidDensity_n( k, q );
177  stack.dFluidMassIncrement_dVolStrainIncrement = dPorosity_dVolStrain * m_fluidDensity( k, q );
178  stack.dFluidMassIncrement_dPressure = dPorosity_dPressure * m_fluidDensity( k, q ) + porosity * m_dFluidDensity( k, q, 0 );
179 }
180 
181 template< typename SUBREGION_TYPE,
182  typename CONSTITUTIVE_TYPE,
183  typename FE_TYPE >
187 assembleMomentumBalanceTerms( real64 const ( &N )[numNodesPerElem],
188  real64 const ( &dNdX )[numNodesPerElem][3],
189  real64 const & detJxW,
190  StackVariables & stack ) const
191 {
192  using namespace PDEUtilities;
193 
194  constexpr FunctionSpace displacementTrialSpace = FE_TYPE::template getFunctionSpace< numDofPerTrialSupportPoint >();
195  constexpr FunctionSpace displacementTestSpace = displacementTrialSpace;
196  constexpr FunctionSpace pressureTrialSpace = FunctionSpace::P0;
197 
198  // Step 1: compute local linear momentum balance residual
199  LinearFormUtilities::compute< displacementTestSpace,
200  DifferentialOperator::SymmetricGradient >
201  (
202  stack.localResidualMomentum,
203  dNdX,
204  stack.totalStress,
205  -detJxW );
206 
207  LinearFormUtilities::compute< displacementTestSpace,
208  DifferentialOperator::Identity >
209  (
210  stack.localResidualMomentum,
211  N,
212  stack.bodyForce,
213  detJxW );
214 
215  // Step 2: compute local linear momentum balance residual derivatives with respect to displacement
216  BilinearFormUtilities::compute< displacementTestSpace,
217  displacementTrialSpace,
218  DifferentialOperator::SymmetricGradient,
219  DifferentialOperator::SymmetricGradient >
220  (
221  stack.dLocalResidualMomentum_dDisplacement,
222  dNdX,
223  stack.stiffness, // fourth-order tensor handled via DiscretizationOps
224  dNdX,
225  -detJxW );
226 
227  BilinearFormUtilities::compute< displacementTestSpace,
228  displacementTrialSpace,
229  DifferentialOperator::Identity,
230  DifferentialOperator::Divergence >
231  (
232  stack.dLocalResidualMomentum_dDisplacement,
233  N,
234  stack.dBodyForce_dVolStrainIncrement,
235  dNdX,
236  detJxW );
237 
238  // Step 3: compute local linear momentum balance residual derivatives with respect to pressure
239  BilinearFormUtilities::compute< displacementTestSpace,
240  pressureTrialSpace,
241  DifferentialOperator::SymmetricGradient,
242  DifferentialOperator::Identity >
243  (
244  stack.dLocalResidualMomentum_dPressure,
245  dNdX,
246  stack.dTotalStress_dPressure,
247  1.0,
248  -detJxW );
249 
250  BilinearFormUtilities::compute< displacementTestSpace,
251  pressureTrialSpace,
252  DifferentialOperator::Identity,
253  DifferentialOperator::Identity >
254  (
255  stack.dLocalResidualMomentum_dPressure,
256  N,
257  stack.dBodyForce_dPressure,
258  1.0,
259  detJxW );
260 }
261 
262 template< typename SUBREGION_TYPE,
263  typename CONSTITUTIVE_TYPE,
264  typename FE_TYPE >
268 assembleElementBasedFlowTerms( real64 const ( &dNdX )[numNodesPerElem][3],
269  real64 const & detJxW,
270  StackVariables & stack ) const
271 {
272  using namespace PDEUtilities;
273 
274  constexpr FunctionSpace displacementTrialSpace = FE_TYPE::template getFunctionSpace< numDofPerTrialSupportPoint >();
275  constexpr FunctionSpace pressureTrialSpace = FunctionSpace::P0;
276  constexpr FunctionSpace pressureTestSpace = pressureTrialSpace;
277 
278  // Step 1: compute local mass balance residual
279  LinearFormUtilities::compute< pressureTestSpace,
280  DifferentialOperator::Identity >
281  (
282  stack.localResidualMass,
283  1.0,
284  stack.fluidMassIncrement,
285  detJxW );
286 
287  // Step 2: compute local mass balance residual derivatives with respect to displacement
288  BilinearFormUtilities::compute< pressureTestSpace,
289  displacementTrialSpace,
290  DifferentialOperator::Identity,
291  DifferentialOperator::Divergence >
292  (
294  1.0,
296  dNdX,
297  detJxW );
298 
299  // Step 3: compute local mass balance residual derivatives with respect to pressure
300  BilinearFormUtilities::compute< pressureTestSpace,
301  pressureTrialSpace,
302  DifferentialOperator::Identity,
303  DifferentialOperator::Identity >
304  (
306  1.0,
308  1.0,
309  detJxW );
310 }
311 
312 template< typename SUBREGION_TYPE,
313  typename CONSTITUTIVE_TYPE,
314  typename FE_TYPE >
319  localIndex const q,
320  StackVariables & stack ) const
321 {
322  // Governing equations (strong form)
323  // ---------------------------------
324  //
325  // divergence( totalStress ) + bodyForce = 0 (quasi-static linear momentum balance)
326  // dFluidMass_dTime + divergence( fluidMassFlux ) = source (fluid phase mass balance)
327  //
328  // with currently the following dependencies on the strainIncrement tensor and pressure
329  //
330  // totalStress = totalStress( strainIncrement, pressure)
331  // bodyForce = bodyForce( strainIncrement, pressure)
332  // fluidMass = fluidMass( strainIncrement, pressure)
333  // fluidMassFlux = fludiMassFlux( pressure)
334  //
335  // Note that the fluidMassFlux will depend on the strainIncrement if a stress-dependent constitutive
336  // model is assumed. A dependency on pressure can also occur in the source term of the mass
337  // balance equation, e.g. if a Peaceman well model is used.
338  //
339  // In this kernel cell-based contributions to Jacobian matrix and residual
340  // vector are computed. The face-based contributions, namely associated with the fluid
341  // mass flux term are computed in a different kernel. The source term in the mass balance
342  // equation is also treated elsewhere.
343  //
344  // Integration in time is performed using a backward Euler scheme (in the mass balance
345  // equation LHS and RHS are multiplied by the timestep).
346  //
347  // Using a weak formulation of the governing equation the following terms are assembled in this kernel
348  //
349  // Rmom = - \int symmetricGradient( \eta ) : totalStress + \int \eta \cdot bodyForce = 0
350  // Rmas = \int \chi ( fluidMass - fluidMass_n) = 0
351  //
352  // dRmom_dVolStrain = - \int_Omega symmetricGradient( \eta ) : dTotalStress_dVolStrain
353  // + \int \eta \cdot dBodyForce_dVolStrain
354  // dRmom_dPressure = - \int_Omega symmetricGradient( \eta ) : dTotalStress_dPressure
355  // + \int \eta \cdot dBodyForce_dPressure
356  // dRmas_dVolStrain = \int \chi dFluidMass_dVolStrain
357  // dRmas_dPressure = \int \chi dFluidMass_dPressure
358  //
359  // with \eta and \chi test basis functions for the displacement and pressure field, respectively.
360  // A continuous interpolation is used for the displacement, with \eta continuous finite element
361  // basis functions. A piecewise-constant approximation is used for the pressure.
362 
363  // Step 1: compute displacement finite element basis functions (N), basis function derivatives (dNdX), and
364  // determinant of the Jacobian transformation matrix times the quadrature weight (detJxW)
365  real64 N[numNodesPerElem]{};
366  real64 dNdX[numNodesPerElem][3]{};
367  FE_TYPE::calcN( q, stack.feStack, N );
368  real64 const detJxW = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.xLocal,
369  stack.feStack, dNdX );
370 
371  // Step 2: compute strain increment
372  LvArray::tensorOps::fill< 6 >( stack.strainIncrement, 0.0 );
373  FE_TYPE::symmetricGradient( dNdX, stack.uhat_local, stack.strainIncrement );
374 
375  // Step 3: compute 1) the total stress, 2) the body force terms, and 3) the fluidMassIncrement
376  // using quantities returned by the PorousSolid constitutive model.
377  // This function also computes the derivatives of these three quantities wrt primary variables
378  smallStrainUpdate( k, q, stack );
379 
380  // Step 4: use the total stress and the body force to increment the local momentum balance residual
381  // This function also fills the local Jacobian rows corresponding to the momentum balance.
382  assembleMomentumBalanceTerms( N, dNdX, detJxW, stack );
383 
384  // Step 5: use the fluid mass increment to increment the local mass balance residual
385  // This function also fills the local Jacobian rows corresponding to the mass balance.
386  assembleElementBasedFlowTerms( dNdX, detJxW, stack );
387 }
388 
389 template< typename SUBREGION_TYPE,
390  typename CONSTITUTIVE_TYPE,
391  typename FE_TYPE >
395 complete( localIndex const k,
396  StackVariables & stack ) const
397 {
398  GEOS_UNUSED_VAR( k );
399  real64 maxForce = 0;
400  localIndex const numSupportPoints =
401  m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
402  integer numDisplacementDofs = numSupportPoints * numDofPerTestSupportPoint;
403 
404  for( int localNode = 0; localNode < numSupportPoints; ++localNode )
405  {
406  for( int dim = 0; dim < numDofPerTestSupportPoint; ++dim )
407  {
408 
409  localIndex const dof = LvArray::integerConversion< localIndex >( stack.localRowDofIndex[numDofPerTestSupportPoint*localNode + dim] - m_dofRankOffset );
410 
411  // we need this check to filter out ghost nodes in the assembly
412  if( dof < 0 || dof >= m_matrix.numRows() )
413  {
414  continue;
415  }
416  m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( dof,
417  stack.localRowDofIndex,
418  stack.dLocalResidualMomentum_dDisplacement[numDofPerTestSupportPoint * localNode + dim],
419  numDisplacementDofs );
420 
421  RAJA::atomicAdd< parallelDeviceAtomic >( &m_rhs[dof], stack.localResidualMomentum[numDofPerTestSupportPoint * localNode + dim] );
422  maxForce = fmax( maxForce, fabs( stack.localResidualMomentum[numDofPerTestSupportPoint * localNode + dim] ) );
423  m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( dof,
424  &stack.localPressureDofIndex,
425  stack.dLocalResidualMomentum_dPressure[numDofPerTestSupportPoint * localNode + dim],
426  1 );
427  }
428  }
429 
430  localIndex const dof = LvArray::integerConversion< localIndex >( stack.localPressureDofIndex - m_dofRankOffset );
431 
432  // we need this check to filter out ghost cells in the assembly
433  if( 0 <= dof && dof < m_matrix.numRows() )
434  {
435  m_matrix.template addToRowBinarySearchUnsorted< serialAtomic >( dof,
436  stack.localRowDofIndex,
438  numDisplacementDofs );
439  m_matrix.template addToRow< serialAtomic >( dof,
440  &stack.localPressureDofIndex,
442  1 );
443  RAJA::atomicAdd< serialAtomic >( &m_rhs[dof], stack.localResidualMass[0] );
444  }
445  return maxForce;
446 }
447 
448 template< typename SUBREGION_TYPE,
449  typename CONSTITUTIVE_TYPE,
450  typename FE_TYPE >
451 template< typename POLICY,
452  typename KERNEL_TYPE >
453 real64
455 kernelLaunch( localIndex const numElems,
456  KERNEL_TYPE const & kernelComponent )
457 {
459 
460  // Define a RAJA reduction variable to get the maximum residual contribution.
461  RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxResidual( 0 );
462 
463  forAll< POLICY >( numElems,
464  [=] GEOS_HOST_DEVICE ( localIndex const k )
465  {
466  typename KERNEL_TYPE::StackVariables stack;
467 
468  kernelComponent.setup( k, stack );
469  for( integer q=0; q<KERNEL_TYPE::numQuadraturePointsPerElem; ++q )
470  {
471  kernelComponent.quadraturePointKernel( k, q, stack );
472  }
473  maxResidual.max( kernelComponent.complete( k, stack ) );
474  } );
475  return maxResidual.get();
476 }
477 
478 
479 } // namespace poromechanicsKernels
480 
481 } // namespace geos
482 
483 #endif // GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_SINGLEPHASEPOROMECHANICS_IMPL_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
#define GEOS_MARK_FUNCTION
Mark function with both Caliper and NVTX if enabled.
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
Implements kernels for solving quasi-static single-phase poromechanics.
GEOS_HOST_DEVICE void smallStrainUpdate(localIndex const k, localIndex const q, StackVariables &stack) const
Helper function to compute 1) the total stress, 2) the body force term, and 3) the fluidMassIncrement...
GEOS_HOST_DEVICE void computeFluidIncrement(localIndex const k, localIndex const q, real64 const &porosity, real64 const &porosity_n, real64 const &dPorosity_dVolStrain, real64 const &dPorosity_dPressure, real64 const &dPorosity_dTemperature, StackVariables &stack) const
Helper function to compute the fluid mass increment and its derivatives wrt primary variables.
static real64 kernelLaunch(localIndex const numElems, KERNEL_TYPE const &kernelComponent)
Kernel Launcher.
GEOS_HOST_DEVICE void computeBodyForce(localIndex const k, localIndex const q, real64 const &porosity, real64 const &dPorosity_dVolStrain, real64 const &dPorosity_dPressure, real64 const &dPorosity_dTemperature, real64 const &dSolidDensity_dPressure, StackVariables &stack) const
Helper function to compute the body force term (\rho g) and its derivatives wrt primary variables.
GEOS_HOST_DEVICE void assembleElementBasedFlowTerms(real64 const (&dNdX)[numNodesPerElem][3], real64 const &detJxW, StackVariables &stack) const
Assemble the local mass balance residual and derivatives using fluid mass/energy increment.
SinglePhasePoromechanics(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, integer const performStressInitialization, string const fluidModelKey)
Constructor.
GEOS_HOST_DEVICE real64 complete(localIndex const k, StackVariables &stack) const
Performs the complete phase for the kernel.
GEOS_HOST_DEVICE void assembleMomentumBalanceTerms(real64 const (&N)[numNodesPerElem], real64 const (&dNdX)[numNodesPerElem][3], real64 const &detJxW, StackVariables &stack) const
Assemble the local linear momentum balance residual and derivatives using total stress and body force...
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
std::string string
String type.
Definition: DataTypes.hpp:90
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
real64 dFluidMassIncrement_dVolStrainIncrement
Derivative of mass accumulation wrt volumetric strain increment.
real64 dLocalResidualMass_dPressure[1][1]
Derivative of mass balance residual wrt pressure.
real64 dLocalResidualMass_dDisplacement[1][Base::StackVariables::numDispDofPerElem]
Derivative of mass balance residual wrt displacement.
real64 dFluidMassIncrement_dPressure
Derivative of mass accumulation wrt pressure.