20 #ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_THERMALSINGLEPHASEPOROMECHANICS_IMPL_HPP_
21 #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_THERMALSINGLEPHASEPOROMECHANICS_IMPL_HPP_
23 #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp"
29 namespace thermalPoromechanicsKernels
32 template<
typename SUBREGION_TYPE,
33 typename CONSTITUTIVE_TYPE,
40 SUBREGION_TYPE
const & elementSubRegion,
41 FE_TYPE
const & finiteElementSpace,
42 CONSTITUTIVE_TYPE & inputConstitutiveType,
48 real64 const (&gravityVector)[3],
49 string const inputFlowDofKey,
50 integer const performStressInitialization,
51 string const fluidModelKey ):
58 inputConstitutiveType,
66 performStressInitialization,
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 >() )
83 template<
typename SUBREGION_TYPE,
84 typename CONSTITUTIVE_TYPE,
92 Base::setup( k, stack );
94 stack.temperature = m_temperature[k];
95 stack.deltaTemperatureFromLastStep = m_temperature[k] - m_temperature_n[k];
98 template<
typename SUBREGION_TYPE,
99 typename CONSTITUTIVE_TYPE,
110 real64 dPorosity_dVolStrain = 0.0;
111 real64 dPorosity_dPressure = 0.0;
112 real64 dPorosity_dTemperature = 0.0;
113 real64 dSolidDensity_dPressure = 0.0;
116 m_constitutiveUpdate.smallStrainUpdatePoromechanics( k, q,
121 stack.deltaTemperatureFromLastStep,
122 stack.strainIncrement,
124 stack.dTotalStress_dPressure,
125 stack.dTotalStress_dTemperature,
127 m_performStressInitialization,
130 dPorosity_dVolStrain,
132 dPorosity_dTemperature,
133 dSolidDensity_dPressure );
136 computeBodyForce( k, q,
138 dPorosity_dVolStrain,
140 dPorosity_dTemperature,
141 dSolidDensity_dPressure,
145 computeFluidIncrement( k, q,
148 dPorosity_dVolStrain,
150 dPorosity_dTemperature,
154 template<
typename SUBREGION_TYPE,
155 typename CONSTITUTIVE_TYPE,
163 real64 const & dPorosity_dVolStrain,
164 real64 const & dPorosity_dPressure,
165 real64 const & dPorosity_dTemperature,
166 real64 const & dSolidDensity_dPressure,
169 Base::computeBodyForce( k, q,
171 dPorosity_dVolStrain,
173 dPorosity_dTemperature,
174 dSolidDensity_dPressure,
177 real64 const dMixtureDens_dTemperature =
178 dPorosity_dTemperature * ( -m_solidDensity( k, q ) + m_fluidDensity( k, q ) )
179 + porosity * m_dFluidDensity_dTemperature( k, q );
181 LvArray::tensorOps::scaledCopy< 3 >( stack.
dBodyForce_dTemperature, m_gravityVector, dMixtureDens_dTemperature );
184 template<
typename SUBREGION_TYPE,
185 typename CONSTITUTIVE_TYPE,
193 real64 const & porosity_n,
194 real64 const & dPorosity_dVolStrain,
195 real64 const & dPorosity_dPressure,
196 real64 const & dPorosity_dTemperature,
200 Base::computeFluidIncrement( k, q,
203 dPorosity_dVolStrain,
205 dPorosity_dTemperature,
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 );
219 + fluidMass * m_dFluidInternalEnergy_dPressure( k, q );
221 + fluidMass * m_dFluidInternalEnergy_dTemperature( k, q );
225 real64 const oneMinusPoro = 1 - porosity;
227 stack.
energyIncrement += oneMinusPoro * m_rockInternalEnergy( k, 0 ) - ( 1 - porosity_n ) * m_rockInternalEnergy_n( k, 0 );
230 stack.
dEnergyIncrement_dTemperature += -dPorosity_dTemperature * m_rockInternalEnergy( k, 0 ) + oneMinusPoro * m_dRockInternalEnergy_dTemperature( k, 0 );
233 template<
typename SUBREGION_TYPE,
234 typename CONSTITUTIVE_TYPE,
240 real64 const ( &dNdX )[numNodesPerElem][3],
244 using namespace PDEUtilities;
246 Base::assembleMomentumBalanceTerms( N, dNdX, detJxW, stack );
248 constexpr FunctionSpace displacementTrialSpace = FE_TYPE::template getFunctionSpace< numDofPerTrialSupportPoint >();
249 constexpr FunctionSpace displacementTestSpace = displacementTrialSpace;
250 constexpr FunctionSpace pressureTrialSpace = FunctionSpace::P0;
254 BilinearFormUtilities::compute< displacementTestSpace,
256 DifferentialOperator::SymmetricGradient,
257 DifferentialOperator::Identity >
261 stack.dTotalStress_dTemperature,
265 BilinearFormUtilities::compute< displacementTestSpace,
267 DifferentialOperator::Identity,
268 DifferentialOperator::Identity >
277 template<
typename SUBREGION_TYPE,
278 typename CONSTITUTIVE_TYPE,
288 using namespace PDEUtilities;
290 Base::assembleElementBasedFlowTerms( dNdX, detJxW, stack );
292 constexpr FunctionSpace displacementTrialSpace = FE_TYPE::template getFunctionSpace< numDofPerTrialSupportPoint >();
293 constexpr FunctionSpace pressureTrialSpace = FunctionSpace::P0;
294 constexpr FunctionSpace pressureTestSpace = pressureTrialSpace;
298 BilinearFormUtilities::compute< pressureTestSpace,
300 DifferentialOperator::Identity,
301 DifferentialOperator::Identity >
311 LinearFormUtilities::compute< pressureTestSpace,
312 DifferentialOperator::Identity >
319 BilinearFormUtilities::compute< pressureTestSpace,
320 displacementTrialSpace,
321 DifferentialOperator::Identity,
322 DifferentialOperator::Divergence >
330 BilinearFormUtilities::compute< pressureTestSpace,
332 DifferentialOperator::Identity,
333 DifferentialOperator::Identity >
341 BilinearFormUtilities::compute< pressureTestSpace,
343 DifferentialOperator::Identity,
344 DifferentialOperator::Identity >
353 template<
typename SUBREGION_TYPE,
354 typename CONSTITUTIVE_TYPE,
361 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,
398 real64 const maxForce = Base::complete( k, stack );
401 m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
402 integer numDisplacementDofs = numSupportPoints * numDofPerTestSupportPoint;
406 for(
int localNode = 0; localNode < numSupportPoints; ++localNode )
408 for(
integer dim = 0; dim < numDofPerTestSupportPoint; ++dim )
410 localIndex const dof = LvArray::integerConversion< localIndex >( stack.localRowDofIndex[numDofPerTestSupportPoint*localNode + dim] - m_dofRankOffset );
413 if( dof < 0 || dof >= m_matrix.numRows() )
418 m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( dof,
427 localIndex const massDof = LvArray::integerConversion< localIndex >( stack.localPressureDofIndex - m_dofRankOffset );
430 if( 0 <= massDof && massDof < m_matrix.numRows() )
432 m_matrix.template addToRow< serialAtomic >( massDof,
443 if( 0 <= energyDof && energyDof < m_matrix.numRows() )
445 m_matrix.template addToRowBinarySearchUnsorted< serialAtomic >( energyDof,
446 stack.localRowDofIndex,
448 numDisplacementDofs );
449 m_matrix.template addToRow< serialAtomic >( energyDof,
450 &stack.localPressureDofIndex,
453 m_matrix.template addToRow< serialAtomic >( energyDof,
464 template<
typename SUBREGION_TYPE,
465 typename CONSTITUTIVE_TYPE,
467 template<
typename POLICY,
468 typename KERNEL_TYPE >
471 KERNEL_TYPE
const & kernelComponent )
476 RAJA::ReduceMax< ReducePolicy< POLICY >,
real64 > maxResidual( 0 );
478 forAll< POLICY >( numElems,
481 typename KERNEL_TYPE::StackVariables stack;
483 kernelComponent.setup( k, stack );
484 for(
integer q=0; q<KERNEL_TYPE::numQuadraturePointsPerElem; ++q )
486 kernelComponent.quadraturePointKernel( k, q, stack );
488 maxResidual.max( kernelComponent.complete( k, stack ) );
490 return maxResidual.get();
#define GEOS_HOST_DEVICE
Marks a host-device function.
#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 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.
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 dLocalResidualEnergy_dDisplacement[1][numDispDofPerElem]
Derivative of energy balance residual wrt displacement.
real64 dEnergyIncrement_dPressure
Derivative of energy accumulation wrt pressure.
real64 dFluidMassIncrement_dTemperature
Derivative of mass accumulation wrt temperature.
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 energyIncrement
Energy accumulation.
real64 dEnergyIncrement_dTemperature
Derivative of energy accumulation wrt temperature.
real64 localResidualEnergy[1]
Energy balance residual.
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 dBodyForce_dTemperature[3]
Derivative of body force wrt temperature.
real64 dLocalResidualEnergy_dTemperature[1][1]
Derivative of energy balance residual wrt temperature.