20 #ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_SINGLEPHASEPOROMECHANICS_IMPL_HPP_
21 #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_SINGLEPHASEPOROMECHANICS_IMPL_HPP_
23 #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp"
24 #include "finiteElement/BilinearFormUtilities.hpp"
25 #include "finiteElement/LinearFormUtilities.hpp"
31 namespace poromechanicsKernels
34 template<
typename SUBREGION_TYPE,
35 typename CONSTITUTIVE_TYPE,
42 SUBREGION_TYPE
const & elementSubRegion,
43 FE_TYPE
const & finiteElementSpace,
44 CONSTITUTIVE_TYPE & inputConstitutiveType,
50 real64 const (&gravityVector)[3],
51 string const inputFlowDofKey,
52 integer const performStressInitialization,
53 string const fluidModelKey ):
60 inputConstitutiveType,
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 )
75 template<
typename SUBREGION_TYPE,
76 typename CONSTITUTIVE_TYPE,
87 real64 dPorosity_dVolStrain = 0.0;
88 real64 dPorosity_dPressure = 0.0;
89 real64 dPorosity_dTemperature = 0.0;
90 real64 dSolidDensity_dPressure = 0.0;
93 m_constitutiveUpdate.smallStrainUpdatePoromechanics( k, q,
98 stack.deltaTemperatureFromLastStep,
99 stack.strainIncrement,
101 stack.dTotalStress_dPressure,
102 stack.dTotalStress_dTemperature,
104 m_performStressInitialization,
107 dPorosity_dVolStrain,
109 dPorosity_dTemperature,
110 dSolidDensity_dPressure );
113 computeBodyForce( k, q,
115 dPorosity_dVolStrain,
117 dPorosity_dTemperature,
118 dSolidDensity_dPressure,
122 computeFluidIncrement( k, q,
125 dPorosity_dVolStrain,
127 dPorosity_dTemperature,
131 template<
typename SUBREGION_TYPE,
132 typename CONSTITUTIVE_TYPE,
140 real64 const & dPorosity_dVolStrain,
141 real64 const & dPorosity_dPressure,
142 real64 const & dPorosity_dTemperature,
143 real64 const & dSolidDensity_dPressure,
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 );
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 );
159 template<
typename SUBREGION_TYPE,
160 typename CONSTITUTIVE_TYPE,
168 real64 const & porosity_n,
169 real64 const & dPorosity_dVolStrain,
170 real64 const & dPorosity_dPressure,
171 real64 const & dPorosity_dTemperature,
176 stack.
fluidMassIncrement = porosity * m_fluidDensity( k, q ) - porosity_n * m_fluidDensity_n( k, q );
181 template<
typename SUBREGION_TYPE,
182 typename CONSTITUTIVE_TYPE,
188 real64 const ( &dNdX )[numNodesPerElem][3],
192 using namespace PDEUtilities;
194 constexpr FunctionSpace displacementTrialSpace = FE_TYPE::template getFunctionSpace< numDofPerTrialSupportPoint >();
195 constexpr FunctionSpace displacementTestSpace = displacementTrialSpace;
196 constexpr FunctionSpace pressureTrialSpace = FunctionSpace::P0;
199 LinearFormUtilities::compute< displacementTestSpace,
200 DifferentialOperator::SymmetricGradient >
202 stack.localResidualMomentum,
207 LinearFormUtilities::compute< displacementTestSpace,
208 DifferentialOperator::Identity >
210 stack.localResidualMomentum,
216 BilinearFormUtilities::compute< displacementTestSpace,
217 displacementTrialSpace,
218 DifferentialOperator::SymmetricGradient,
219 DifferentialOperator::SymmetricGradient >
221 stack.dLocalResidualMomentum_dDisplacement,
227 BilinearFormUtilities::compute< displacementTestSpace,
228 displacementTrialSpace,
229 DifferentialOperator::Identity,
230 DifferentialOperator::Divergence >
232 stack.dLocalResidualMomentum_dDisplacement,
234 stack.dBodyForce_dVolStrainIncrement,
239 BilinearFormUtilities::compute< displacementTestSpace,
241 DifferentialOperator::SymmetricGradient,
242 DifferentialOperator::Identity >
244 stack.dLocalResidualMomentum_dPressure,
246 stack.dTotalStress_dPressure,
250 BilinearFormUtilities::compute< displacementTestSpace,
252 DifferentialOperator::Identity,
253 DifferentialOperator::Identity >
255 stack.dLocalResidualMomentum_dPressure,
257 stack.dBodyForce_dPressure,
262 template<
typename SUBREGION_TYPE,
263 typename CONSTITUTIVE_TYPE,
272 using namespace PDEUtilities;
274 constexpr FunctionSpace displacementTrialSpace = FE_TYPE::template getFunctionSpace< numDofPerTrialSupportPoint >();
275 constexpr FunctionSpace pressureTrialSpace = FunctionSpace::P0;
276 constexpr FunctionSpace pressureTestSpace = pressureTrialSpace;
279 LinearFormUtilities::compute< pressureTestSpace,
280 DifferentialOperator::Identity >
288 BilinearFormUtilities::compute< pressureTestSpace,
289 displacementTrialSpace,
290 DifferentialOperator::Identity,
291 DifferentialOperator::Divergence >
300 BilinearFormUtilities::compute< pressureTestSpace,
302 DifferentialOperator::Identity,
303 DifferentialOperator::Identity >
312 template<
typename SUBREGION_TYPE,
313 typename CONSTITUTIVE_TYPE,
320 StackVariables & stack )
const
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 );
372 LvArray::tensorOps::fill< 6 >( stack.strainIncrement, 0.0 );
373 FE_TYPE::symmetricGradient( dNdX, stack.uhat_local, stack.strainIncrement );
378 smallStrainUpdate( k, q, stack );
382 assembleMomentumBalanceTerms( N, dNdX, detJxW, stack );
386 assembleElementBasedFlowTerms( dNdX, detJxW, stack );
389 template<
typename SUBREGION_TYPE,
390 typename CONSTITUTIVE_TYPE,
401 m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
402 integer numDisplacementDofs = numSupportPoints * numDofPerTestSupportPoint;
404 for(
int localNode = 0; localNode < numSupportPoints; ++localNode )
406 for(
int dim = 0; dim < numDofPerTestSupportPoint; ++dim )
409 localIndex const dof = LvArray::integerConversion< localIndex >( stack.localRowDofIndex[numDofPerTestSupportPoint*localNode + dim] - m_dofRankOffset );
412 if( dof < 0 || dof >= m_matrix.numRows() )
416 m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( dof,
417 stack.localRowDofIndex,
418 stack.dLocalResidualMomentum_dDisplacement[numDofPerTestSupportPoint * localNode + dim],
419 numDisplacementDofs );
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],
430 localIndex const dof = LvArray::integerConversion< localIndex >( stack.localPressureDofIndex - m_dofRankOffset );
433 if( 0 <= dof && dof < m_matrix.numRows() )
435 m_matrix.template addToRowBinarySearchUnsorted< serialAtomic >( dof,
436 stack.localRowDofIndex,
438 numDisplacementDofs );
439 m_matrix.template addToRow< serialAtomic >( dof,
440 &stack.localPressureDofIndex,
448 template<
typename SUBREGION_TYPE,
449 typename CONSTITUTIVE_TYPE,
451 template<
typename POLICY,
452 typename KERNEL_TYPE >
456 KERNEL_TYPE
const & kernelComponent )
461 RAJA::ReduceMax< ReducePolicy< POLICY >,
real64 > maxResidual( 0 );
463 forAll< POLICY >( numElems,
466 typename KERNEL_TYPE::StackVariables stack;
468 kernelComponent.setup( k, stack );
469 for(
integer q=0; q<KERNEL_TYPE::numQuadraturePointsPerElem; ++q )
471 kernelComponent.quadraturePointKernel( k, q, stack );
473 maxResidual.max( kernelComponent.complete( k, stack ) );
475 return maxResidual.get();
#define GEOS_HOST_DEVICE
Marks a host-device function.
#define GEOS_UNUSED_VAR(...)
Mark an unused variable and silence compiler warnings.
#define GEOS_FORCE_INLINE
Marks a function or lambda for inlining.
#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.
The FaceManager class provides an interface to ObjectManagerBase in order to manage face data.
The NodeManager class provides an interface to ObjectManagerBase in order to manage node data.
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.
LvArray::CRSMatrixView< T, COL_INDEX, localIndex const, LvArray::ChaiBuffer > CRSMatrixView
Alias for CRS Matrix View.
GEOS_GLOBALINDEX_TYPE globalIndex
Global index type (for indexing objects across MPI partitions).
std::string string
String type.
double real64
64-bit floating point type.
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
std::int32_t integer
Signed integer type.
Kernel variables allocated on the stack.
real64 dFluidMassIncrement_dVolStrainIncrement
Derivative of mass accumulation wrt volumetric strain increment.
real64 fluidMassIncrement
Mass accumulation.
real64 localResidualMass[1]
Mass balance residual.
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.