GEOS
ThermalSinglePhasePoromechanics_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_THERMALSINGLEPHASEPOROMECHANICS_IMPL_HPP_
21 #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_THERMALSINGLEPHASEPOROMECHANICS_IMPL_HPP_
22 
23 #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp"
25 
26 namespace geos
27 {
28 
29 namespace thermalPoromechanicsKernels
30 {
31 
32 template< typename SUBREGION_TYPE,
33  typename CONSTITUTIVE_TYPE,
34  typename FE_TYPE >
36 ThermalSinglePhasePoromechanics( NodeManager const & nodeManager,
37  EdgeManager const & edgeManager,
38  FaceManager const & faceManager,
39  localIndex const targetRegionIndex,
40  SUBREGION_TYPE const & elementSubRegion,
41  FE_TYPE const & finiteElementSpace,
42  CONSTITUTIVE_TYPE & inputConstitutiveType,
43  arrayView1d< globalIndex const > const inputDispDofNumber,
44  globalIndex const rankOffset,
46  arrayView1d< real64 > const inputRhs,
47  real64 const inputDt,
48  real64 const (&gravityVector)[3],
49  string const inputFlowDofKey,
50  integer const performStressInitialization,
51  string const fluidModelKey ):
52  Base( nodeManager,
53  edgeManager,
54  faceManager,
55  targetRegionIndex,
56  elementSubRegion,
57  finiteElementSpace,
58  inputConstitutiveType,
59  inputDispDofNumber,
60  rankOffset,
61  inputMatrix,
62  inputRhs,
63  inputDt,
64  gravityVector,
65  inputFlowDofKey,
66  performStressInitialization,
67  fluidModelKey ),
68  m_dFluidDensity_dTemperature( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >(
69  fluidModelKey ) ).dDensity_dTemperature() ),
70  m_fluidInternalEnergy_n( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).internalEnergy_n() ),
71  m_fluidInternalEnergy( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).internalEnergy() ),
72  m_dFluidInternalEnergy_dPressure( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >(
73  fluidModelKey ) ).dInternalEnergy_dPressure() ),
74  m_dFluidInternalEnergy_dTemperature( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >(
75  fluidModelKey ) ).dInternalEnergy_dTemperature() ),
76  m_rockInternalEnergy_n( inputConstitutiveType.getInternalEnergy_n() ),
77  m_rockInternalEnergy( inputConstitutiveType.getInternalEnergy() ),
78  m_dRockInternalEnergy_dTemperature( inputConstitutiveType.getDinternalEnergy_dTemperature() ),
79  m_temperature_n( elementSubRegion.template getField< fields::flow::temperature_n >() ),
80  m_temperature( elementSubRegion.template getField< fields::flow::temperature >() )
81 {}
82 
83 template< typename SUBREGION_TYPE,
84  typename CONSTITUTIVE_TYPE,
85  typename FE_TYPE >
89 setup( localIndex const k,
90  StackVariables & stack ) const
91 {
92  Base::setup( k, stack );
93  stack.localTemperatureDofIndex = m_flowDofNumber[k]+1;
94  stack.temperature = m_temperature[k];
95  stack.deltaTemperatureFromLastStep = m_temperature[k] - m_temperature_n[k];
96 }
97 
98 template< typename SUBREGION_TYPE,
99  typename CONSTITUTIVE_TYPE,
100  typename FE_TYPE >
105  localIndex const q,
106  StackVariables & stack ) const
107 {
108  real64 porosity = 0.0;
109  real64 porosity_n = 0.0;
110  real64 dPorosity_dVolStrain = 0.0;
111  real64 dPorosity_dPressure = 0.0;
112  real64 dPorosity_dTemperature = 0.0;
113  real64 dSolidDensity_dPressure = 0.0;
114 
115  // Step 1: call the constitutive model to evaluate the total stress and compute porosity
116  m_constitutiveUpdate.smallStrainUpdatePoromechanics( k, q,
117  m_dt,
118  m_pressure[k],
119  m_pressure_n[k],
120  stack.temperature,
121  stack.deltaTemperatureFromLastStep,
122  stack.strainIncrement,
123  stack.totalStress,
124  stack.dTotalStress_dPressure,
125  stack.dTotalStress_dTemperature,
126  stack.stiffness,
127  m_performStressInitialization,
128  porosity,
129  porosity_n,
130  dPorosity_dVolStrain,
131  dPorosity_dPressure,
132  dPorosity_dTemperature,
133  dSolidDensity_dPressure );
134 
135  // Step 2: compute the body force
136  computeBodyForce( k, q,
137  porosity,
138  dPorosity_dVolStrain,
139  dPorosity_dPressure,
140  dPorosity_dTemperature,
141  dSolidDensity_dPressure,
142  stack );
143 
144  // Step 3: compute fluid mass increment
145  computeFluidIncrement( k, q,
146  porosity,
147  porosity_n,
148  dPorosity_dVolStrain,
149  dPorosity_dPressure,
150  dPorosity_dTemperature,
151  stack );
152 }
153 
154 template< typename SUBREGION_TYPE,
155  typename CONSTITUTIVE_TYPE,
156  typename FE_TYPE >
161  localIndex const q,
162  real64 const & porosity,
163  real64 const & dPorosity_dVolStrain,
164  real64 const & dPorosity_dPressure,
165  real64 const & dPorosity_dTemperature,
166  real64 const & dSolidDensity_dPressure,
167  StackVariables & stack ) const
168 {
169  Base::computeBodyForce( k, q,
170  porosity,
171  dPorosity_dVolStrain,
172  dPorosity_dPressure,
173  dPorosity_dTemperature,
174  dSolidDensity_dPressure,
175  stack );
176 
177  real64 const dMixtureDens_dTemperature =
178  dPorosity_dTemperature * ( -m_solidDensity( k, q ) + m_fluidDensity( k, q ) )
179  + porosity * m_dFluidDensity_dTemperature( k, q );
180 
181  LvArray::tensorOps::scaledCopy< 3 >( stack.dBodyForce_dTemperature, m_gravityVector, dMixtureDens_dTemperature );
182 }
183 
184 template< typename SUBREGION_TYPE,
185  typename CONSTITUTIVE_TYPE,
186  typename FE_TYPE >
191  localIndex const q,
192  real64 const & porosity,
193  real64 const & porosity_n,
194  real64 const & dPorosity_dVolStrain,
195  real64 const & dPorosity_dPressure,
196  real64 const & dPorosity_dTemperature,
197  StackVariables & stack ) const
198 {
199  // Step 1: compute fluid mass increment and its derivatives wrt vol strain and pressure
200  Base::computeFluidIncrement( k, q,
201  porosity,
202  porosity_n,
203  dPorosity_dVolStrain,
204  dPorosity_dPressure,
205  dPorosity_dTemperature,
206  stack );
207 
208  // Step 2: compute derivative of fluid mass increment wrt temperature
209  stack.dFluidMassIncrement_dTemperature = dPorosity_dTemperature * m_fluidDensity( k, q ) + porosity * m_dFluidDensity_dTemperature( k, q );
210 
211  // Step 3: compute fluid energy increment and its derivatives wrt vol strain, pressure, and temperature
212  real64 const fluidMass = porosity * m_fluidDensity( k, q );
213  real64 const fluidEnergy = fluidMass * m_fluidInternalEnergy( k, q );
214  real64 const fluidEnergy_n = porosity_n * m_fluidDensity_n( k, q ) * m_fluidInternalEnergy_n( k, q );
215  stack.energyIncrement = fluidEnergy - fluidEnergy_n;
216 
217  stack.dEnergyIncrement_dVolStrainIncrement = stack.dFluidMassIncrement_dVolStrainIncrement * m_fluidInternalEnergy( k, q );
218  stack.dEnergyIncrement_dPressure = stack.dFluidMassIncrement_dPressure * m_fluidInternalEnergy( k, q )
219  + fluidMass * m_dFluidInternalEnergy_dPressure( k, q );
220  stack.dEnergyIncrement_dTemperature = stack.dFluidMassIncrement_dTemperature * m_fluidInternalEnergy( k, q )
221  + fluidMass * m_dFluidInternalEnergy_dTemperature( k, q );
222 
223 
224  // Step 4: assemble the solid part of the accumulation term
225  real64 const oneMinusPoro = 1 - porosity;
226 
227  stack.energyIncrement += oneMinusPoro * m_rockInternalEnergy( k, 0 ) - ( 1 - porosity_n ) * m_rockInternalEnergy_n( k, 0 );
228  stack.dEnergyIncrement_dVolStrainIncrement += -dPorosity_dVolStrain * m_rockInternalEnergy( k, 0 );
229  stack.dEnergyIncrement_dPressure += -dPorosity_dPressure * m_rockInternalEnergy( k, 0 );
230  stack.dEnergyIncrement_dTemperature += -dPorosity_dTemperature * m_rockInternalEnergy( k, 0 ) + oneMinusPoro * m_dRockInternalEnergy_dTemperature( k, 0 );
231 }
232 
233 template< typename SUBREGION_TYPE,
234  typename CONSTITUTIVE_TYPE,
235  typename FE_TYPE >
239 assembleMomentumBalanceTerms( real64 const ( &N )[numNodesPerElem],
240  real64 const ( &dNdX )[numNodesPerElem][3],
241  real64 const & detJxW,
242  StackVariables & stack ) const
243 {
244  using namespace PDEUtilities;
245 
246  Base::assembleMomentumBalanceTerms( N, dNdX, detJxW, stack );
247 
248  constexpr FunctionSpace displacementTrialSpace = FE_TYPE::template getFunctionSpace< numDofPerTrialSupportPoint >();
249  constexpr FunctionSpace displacementTestSpace = displacementTrialSpace;
250  constexpr FunctionSpace pressureTrialSpace = FunctionSpace::P0;
251 
252  // compute local linear momentum balance residual derivatives with respect to temperature
253 
254  BilinearFormUtilities::compute< displacementTestSpace,
255  pressureTrialSpace,
256  DifferentialOperator::SymmetricGradient,
257  DifferentialOperator::Identity >
258  (
260  dNdX,
261  stack.dTotalStress_dTemperature,
262  1.0,
263  -detJxW );
264 
265  BilinearFormUtilities::compute< displacementTestSpace,
266  pressureTrialSpace,
267  DifferentialOperator::Identity,
268  DifferentialOperator::Identity >
269  (
271  N,
273  1.0,
274  detJxW );
275 }
276 
277 template< typename SUBREGION_TYPE,
278  typename CONSTITUTIVE_TYPE,
279  typename FE_TYPE >
283 assembleElementBasedFlowTerms( real64 const ( &dNdX )[numNodesPerElem][3],
284  real64 const & detJxW,
285  StackVariables & stack ) const
286 {
287 
288  using namespace PDEUtilities;
289 
290  Base::assembleElementBasedFlowTerms( dNdX, detJxW, stack );
291 
292  constexpr FunctionSpace displacementTrialSpace = FE_TYPE::template getFunctionSpace< numDofPerTrialSupportPoint >();
293  constexpr FunctionSpace pressureTrialSpace = FunctionSpace::P0;
294  constexpr FunctionSpace pressureTestSpace = pressureTrialSpace;
295 
296  // Step 1: compute local mass balance residual derivatives with respect to temperature
297 
298  BilinearFormUtilities::compute< pressureTestSpace,
299  pressureTrialSpace,
300  DifferentialOperator::Identity,
301  DifferentialOperator::Identity >
302  (
304  1.0,
306  1.0,
307  detJxW );
308 
309  // Step 2: compute local energy balance residual and its derivatives
310 
311  LinearFormUtilities::compute< pressureTestSpace,
312  DifferentialOperator::Identity >
313  (
314  stack.localResidualEnergy,
315  1.0,
316  stack.energyIncrement,
317  detJxW );
318 
319  BilinearFormUtilities::compute< pressureTestSpace,
320  displacementTrialSpace,
321  DifferentialOperator::Identity,
322  DifferentialOperator::Divergence >
323  (
325  1.0,
327  dNdX,
328  detJxW );
329 
330  BilinearFormUtilities::compute< pressureTestSpace,
331  pressureTrialSpace,
332  DifferentialOperator::Identity,
333  DifferentialOperator::Identity >
334  (
336  1.0,
338  1.0,
339  detJxW );
340 
341  BilinearFormUtilities::compute< pressureTestSpace,
342  pressureTrialSpace,
343  DifferentialOperator::Identity,
344  DifferentialOperator::Identity >
345  (
347  1.0,
349  1.0,
350  detJxW );
351 }
352 
353 template< typename SUBREGION_TYPE,
354  typename CONSTITUTIVE_TYPE,
355  typename FE_TYPE >
360  localIndex const q,
361  StackVariables & stack ) const
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  real64 const maxForce = Base::complete( k, stack );
399 
400  localIndex const numSupportPoints =
401  m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
402  integer numDisplacementDofs = numSupportPoints * numDofPerTestSupportPoint;
403 
404  // Step 1: assemble the derivatives of linear momentum balance wrt temperature into the global matrix
405 
406  for( int localNode = 0; localNode < numSupportPoints; ++localNode )
407  {
408  for( integer dim = 0; dim < numDofPerTestSupportPoint; ++dim )
409  {
410  localIndex const dof = LvArray::integerConversion< localIndex >( stack.localRowDofIndex[numDofPerTestSupportPoint*localNode + dim] - m_dofRankOffset );
411 
412  // we need this check to filter out ghost nodes in the assembly
413  if( dof < 0 || dof >= m_matrix.numRows() )
414  {
415  continue;
416  }
417 
418  m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( dof,
420  stack.dLocalResidualMomentum_dTemperature[numDofPerTestSupportPoint * localNode + dim],
421  1 );
422  }
423  }
424 
425  // Step 2: assemble the derivatives of mass balance residual wrt temperature into the global matrix
426 
427  localIndex const massDof = LvArray::integerConversion< localIndex >( stack.localPressureDofIndex - m_dofRankOffset );
428 
429  // we need this check to filter out ghost cells in the assembly
430  if( 0 <= massDof && massDof < m_matrix.numRows() )
431  {
432  m_matrix.template addToRow< serialAtomic >( massDof,
435  1 );
436  }
437 
438  // Step 3: assemble the energy balance and its derivatives into the global matrix
439 
440  localIndex const energyDof = LvArray::integerConversion< localIndex >( stack.localTemperatureDofIndex - m_dofRankOffset );
441 
442  // we need this check to filter out ghost cells in the assembly
443  if( 0 <= energyDof && energyDof < m_matrix.numRows() )
444  {
445  m_matrix.template addToRowBinarySearchUnsorted< serialAtomic >( energyDof,
446  stack.localRowDofIndex,
448  numDisplacementDofs );
449  m_matrix.template addToRow< serialAtomic >( energyDof,
450  &stack.localPressureDofIndex,
452  1 );
453  m_matrix.template addToRow< serialAtomic >( energyDof,
456  1 );
457 
458  RAJA::atomicAdd< serialAtomic >( &m_rhs[energyDof], stack.localResidualEnergy[0] );
459  }
460 
461  return maxForce;
462 }
463 
464 template< typename SUBREGION_TYPE,
465  typename CONSTITUTIVE_TYPE,
466  typename FE_TYPE >
467 template< typename POLICY,
468  typename KERNEL_TYPE >
470 kernelLaunch( localIndex const numElems,
471  KERNEL_TYPE const & kernelComponent )
472 {
474 
475  // Define a RAJA reduction variable to get the maximum residual contribution.
476  RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxResidual( 0 );
477 
478  forAll< POLICY >( numElems,
479  [=] GEOS_HOST_DEVICE ( localIndex const k )
480  {
481  typename KERNEL_TYPE::StackVariables stack;
482 
483  kernelComponent.setup( k, stack );
484  for( integer q=0; q<KERNEL_TYPE::numQuadraturePointsPerElem; ++q )
485  {
486  kernelComponent.quadraturePointKernel( k, q, stack );
487  }
488  maxResidual.max( kernelComponent.complete( k, stack ) );
489  } );
490  return maxResidual.get();
491 }
492 
493 
494 } // namespace thermalPoromechanicsKernels
495 
496 } // namespace geos
497 
498 #endif // GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_THERMALSINGLEPHASEPOROMECHANICS_HPP_
#define GEOS_HOST_DEVICE
Marks a host-device function.
Definition: GeosxMacros.hpp:49
#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 thermal single-phase poromechanics.
GEOS_HOST_DEVICE void assembleElementBasedFlowTerms(real64 const (&dNdX)[numNodesPerElem][3], real64 const &detJxW, StackVariables &stack) const
Assemble the local mass/energy balance residual and derivatives using fluid mass/energy increment.
static real64 kernelLaunch(localIndex const numElems, KERNEL_TYPE const &kernelComponent)
Kernel Launcher.
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/energy increment and its derivatives wrt primary variables.
GEOS_HOST_DEVICE void setup(localIndex const k, StackVariables &stack) const
Copy global values from primary field to a local stack array.
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 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...
ThermalSinglePhasePoromechanics(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 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 fluidMass/EnergyIn...
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
Definition: DataTypes.hpp:180
LvArray::CRSMatrixView< T, COL_INDEX, localIndex const, LvArray::ChaiBuffer > CRSMatrixView
Alias for CRS Matrix View.
Definition: DataTypes.hpp:310
GEOS_GLOBALINDEX_TYPE globalIndex
Global index type (for indexing objects across MPI partitions).
Definition: DataTypes.hpp:88
std::string string
String type.
Definition: DataTypes.hpp:91
double real64
64-bit floating point type.
Definition: DataTypes.hpp:99
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Definition: DataTypes.hpp:85
std::int32_t integer
Signed integer type.
Definition: DataTypes.hpp:82
real64 dLocalResidualEnergy_dDisplacement[1][numDispDofPerElem]
Derivative of energy balance residual wrt displacement.
globalIndex localTemperatureDofIndex
C-array storage for the element local row degrees of freedom.
real64 dEnergyIncrement_dVolStrainIncrement
Derivative of energy accumulation wrt volumetric strain increment.
real64 dLocalResidualMass_dTemperature[1][1]
Derivative of mass balance residual wrt pressure.
real64 dLocalResidualEnergy_dPressure[1][1]
Derivative of energy balance residual wrt pressure.
real64 dLocalResidualMomentum_dTemperature[numDispDofPerElem][1]
Derivative of linear momentum balance residual wrt temperature.
real64 dLocalResidualEnergy_dTemperature[1][1]
Derivative of energy balance residual wrt temperature.