20 #ifndef GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_FIXEDSTRESSTHERMOPOROMECHANICS_IMPL_HPP_
21 #define GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_FIXEDSTRESSTHERMOPOROMECHANICS_IMPL_HPP_
31 namespace solidMechanicsLagrangianFEMKernels
34 template<
typename SUBREGION_TYPE,
35 typename CONSTITUTIVE_TYPE,
43 SUBREGION_TYPE
const & elementSubRegion,
44 FE_TYPE
const & finiteElementSpace,
45 CONSTITUTIVE_TYPE & inputConstitutiveType,
51 real64 const (&inputGravityVector)[3] ):
58 inputConstitutiveType,
64 m_X( nodeManager.referencePosition()),
65 m_disp( nodeManager.getField< fields::solidMechanics::totalDisplacement >() ),
66 m_uhat( nodeManager.getField< fields::solidMechanics::incrementalDisplacement >() ),
67 m_gravityVector{ inputGravityVector[0], inputGravityVector[1], inputGravityVector[2] },
68 m_bulkDensity( elementSubRegion.template getField< fields::poromechanics::bulkDensity >() ),
69 m_pressure( elementSubRegion.template getField< fields::flow::pressure >() ),
70 m_pressure_n( elementSubRegion.template getField< fields::flow::pressure_n >() ),
71 m_initialTemperature( elementSubRegion.template getField< fields::flow::initialTemperature >() ),
72 m_temperature( elementSubRegion.template getField< fields::flow::temperature >() ),
73 m_temperature_n( elementSubRegion.template getField< fields::flow::temperature_n >() )
76 template<
typename SUBREGION_TYPE,
77 typename CONSTITUTIVE_TYPE,
85 m_finiteElementSpace.template setup< FE_TYPE >( k, m_meshData, stack.feStack );
87 m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
88 stack.numRows = 3 * numSupportPoints;
89 stack.numCols = stack.numRows;
91 for(
localIndex a = 0; a < numSupportPoints; ++a )
93 localIndex const localNodeIndex = m_elemsToNodes( k, a );
95 for(
int i = 0; i < 3; ++i )
97 #if defined(CALC_FEM_SHAPE_IN_KERNEL)
98 stack.
xLocal[ a ][ i ] = m_X[ localNodeIndex ][ i ];
100 stack.
u_local[ a ][i] = m_disp[ localNodeIndex ][i];
101 stack.
uhat_local[ a ][i] = m_uhat[ localNodeIndex ][i];
102 stack.localRowDofIndex[a*3+i] = m_dofNumber[localNodeIndex]+i;
103 stack.localColDofIndex[a*3+i] = m_dofNumber[localNodeIndex]+i;
109 real64 const stabilizationScaling = computeStabilizationScaling( k );
110 m_finiteElementSpace.template addGradGradStabilizationMatrix
111 < FE_TYPE, numDofPerTrialSupportPoint,
true >( stack.feStack,
113 -stabilizationScaling );
116 template<
typename SUBREGION_TYPE,
117 typename CONSTITUTIVE_TYPE,
126 real64 dNdX[ numNodesPerElem ][ 3 ];
127 real64 const detJxW = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.
xLocal,
128 stack.feStack, dNdX );
130 real64 strainInc[6] = {0};
131 real64 totalStress[6] = {0};
133 typename CONSTITUTIVE_TYPE::KernelWrapper::DiscretizationOps stiffness;
135 FE_TYPE::symmetricGradient( dNdX, stack.
uhat_local, strainInc );
140 m_constitutiveUpdate.smallStrainUpdatePoromechanicsFixedStress( k, q,
146 m_initialTemperature[k],
153 totalStress[i] *= -detJxW;
158 real64 const gravityForce[3] = { m_gravityVector[0] * m_bulkDensity( k, q )* detJxW,
159 m_gravityVector[1] * m_bulkDensity( k, q )* detJxW,
160 m_gravityVector[2] * m_bulkDensity( k, q )* detJxW };
162 real64 N[numNodesPerElem];
163 FE_TYPE::calcN( q, stack.feStack, N );
164 FE_TYPE::plusGradNajAijPlusNaFi( dNdX,
168 reinterpret_cast< real64 (&)[numNodesPerElem][3]
>(stack.localResidual) );
169 real64 const stabilizationScaling = computeStabilizationScaling( k );
170 m_finiteElementSpace.template
171 addEvaluatedGradGradStabilizationVector< FE_TYPE,
172 numDofPerTrialSupportPoint >( stack.feStack,
174 reinterpret_cast< real64 (&)[numNodesPerElem][3]
>(stack.localResidual),
175 -stabilizationScaling );
176 stiffness.template upperBTDB< numNodesPerElem >( dNdX, -detJxW, stack.localJacobian );
179 template<
typename SUBREGION_TYPE,
180 typename CONSTITUTIVE_TYPE,
192 CONSTITUTIVE_TYPE::KernelWrapper::DiscretizationOps::template fillLowerBTDB< numNodesPerElem >( stack.localJacobian );
194 m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
195 for(
int localNode = 0; localNode < numSupportPoints; ++localNode )
197 for(
int dim = 0; dim < numDofPerTestSupportPoint; ++dim )
200 LvArray::integerConversion< localIndex >( stack.localRowDofIndex[ numDofPerTestSupportPoint * localNode + dim ] - m_dofRankOffset );
201 if( dof < 0 || dof >= m_matrix.numRows() )
203 m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( dof,
204 stack.localRowDofIndex,
205 stack.localJacobian[ numDofPerTestSupportPoint * localNode + dim ],
208 RAJA::atomicAdd< parallelDeviceAtomic >( &m_rhs[ dof ], stack.localResidual[ numDofPerTestSupportPoint * localNode + dim ] );
209 maxForce = fmax( maxForce, fabs( stack.localResidual[ numDofPerTestSupportPoint * localNode + dim ] ) );
215 template<
typename SUBREGION_TYPE,
216 typename CONSTITUTIVE_TYPE,
218 template<
typename POLICY,
219 typename KERNEL_TYPE >
222 KERNEL_TYPE
const & kernelComponent )
224 return Base::template kernelLaunch< POLICY, KERNEL_TYPE >( numElems, kernelComponent );
#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.
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.
static real64 kernelLaunch(localIndex const numElems, KERNEL_TYPE const &kernelComponent)
Kernel Launcher.
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 GEOS_FORCE_INLINE real64 complete(localIndex const k, StackVariables &stack) const
Performs the complete phase for the kernel.
FixedStressThermoPoromechanics(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 inputDofNumber, globalIndex const rankOffset, CRSMatrixView< real64, globalIndex const > const inputMatrix, arrayView1d< real64 > const inputRhs, real64 const inputDt, real64 const (&inputGravityVector)[3])
Constructor.
GEOS_HOST_DEVICE void quadraturePointKernel(localIndex const k, localIndex const q, StackVariables &stack) const
Performs a state update at a quadrature point.
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).
double real64
64-bit floating point type.
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Kernel variables allocated on the stack.
real64 u_local[numNodesPerElem][numDofPerTrialSupportPoint]
Stack storage for the element local nodal displacement.
real64 uhat_local[numNodesPerElem][numDofPerTrialSupportPoint]
Stack storage for the element local nodal incremental displacement.