19 #ifndef GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_FIXEDSTRESSTHERMOPOROMECHANICS_IMPL_HPP_
20 #define GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_FIXEDSTRESSTHERMOPOROMECHANICS_IMPL_HPP_
23 #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp"
24 #include "physicsSolvers/multiphysics/PoromechanicsFields.hpp"
30 namespace solidMechanicsLagrangianFEMKernels
33 template<
typename SUBREGION_TYPE,
34 typename CONSTITUTIVE_TYPE,
42 SUBREGION_TYPE
const & elementSubRegion,
43 FE_TYPE
const & finiteElementSpace,
44 CONSTITUTIVE_TYPE & inputConstitutiveType,
50 real64 const (&inputGravityVector)[3] ):
57 inputConstitutiveType,
63 m_X( nodeManager.referencePosition()),
64 m_disp( nodeManager.getField< fields::solidMechanics::totalDisplacement >() ),
65 m_uhat( nodeManager.getField< fields::solidMechanics::incrementalDisplacement >() ),
66 m_gravityVector{ inputGravityVector[0], inputGravityVector[1], inputGravityVector[2] },
67 m_bulkDensity( elementSubRegion.template getField< fields::poromechanics::bulkDensity >() ),
68 m_pressure( elementSubRegion.template getField< fields::flow::pressure >() ),
69 m_pressure_n( elementSubRegion.template getField< fields::flow::pressure_n >() ),
70 m_initialTemperature( elementSubRegion.template getField< fields::flow::initialTemperature >() ),
71 m_temperature( elementSubRegion.template getField< fields::flow::temperature >() ),
72 m_temperature_n( elementSubRegion.template getField< fields::flow::temperature_n >() )
75 template<
typename SUBREGION_TYPE,
76 typename CONSTITUTIVE_TYPE,
84 m_finiteElementSpace.template setup< FE_TYPE >( k, m_meshData, stack.feStack );
86 m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
87 stack.numRows = 3 * numSupportPoints;
88 stack.numCols = stack.numRows;
90 for(
localIndex a = 0; a < numSupportPoints; ++a )
92 localIndex const localNodeIndex = m_elemsToNodes( k, a );
94 for(
int i = 0; i < 3; ++i )
96 #if defined(CALC_FEM_SHAPE_IN_KERNEL)
97 stack.
xLocal[ a ][ i ] = m_X[ localNodeIndex ][ i ];
99 stack.
u_local[ a ][i] = m_disp[ localNodeIndex ][i];
100 stack.
uhat_local[ a ][i] = m_uhat[ localNodeIndex ][i];
101 stack.localRowDofIndex[a*3+i] = m_dofNumber[localNodeIndex]+i;
102 stack.localColDofIndex[a*3+i] = m_dofNumber[localNodeIndex]+i;
108 real64 const stabilizationScaling = computeStabilizationScaling( k );
109 m_finiteElementSpace.template addGradGradStabilizationMatrix
110 < FE_TYPE, numDofPerTrialSupportPoint,
true >( stack.feStack,
112 -stabilizationScaling );
115 template<
typename SUBREGION_TYPE,
116 typename CONSTITUTIVE_TYPE,
125 real64 dNdX[ numNodesPerElem ][ 3 ];
126 real64 const detJxW = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.
xLocal,
127 stack.feStack, dNdX );
129 real64 strainInc[6] = {0};
130 real64 totalStress[6] = {0};
132 typename CONSTITUTIVE_TYPE::KernelWrapper::DiscretizationOps stiffness;
134 FE_TYPE::symmetricGradient( dNdX, stack.
uhat_local, strainInc );
139 m_constitutiveUpdate.smallStrainUpdatePoromechanicsFixedStress( k, q,
151 totalStress[i] *= -detJxW;
156 real64 const gravityForce[3] = { m_gravityVector[0] * m_bulkDensity( k, q )* detJxW,
157 m_gravityVector[1] * m_bulkDensity( k, q )* detJxW,
158 m_gravityVector[2] * m_bulkDensity( k, q )* detJxW };
160 real64 N[numNodesPerElem];
161 FE_TYPE::calcN( q, stack.feStack, N );
162 FE_TYPE::plusGradNajAijPlusNaFi( dNdX,
166 reinterpret_cast< real64 (&)[numNodesPerElem][3]
>(stack.localResidual) );
167 real64 const stabilizationScaling = computeStabilizationScaling( k );
168 m_finiteElementSpace.template
169 addEvaluatedGradGradStabilizationVector< FE_TYPE,
170 numDofPerTrialSupportPoint >( stack.feStack,
172 reinterpret_cast< real64 (&)[numNodesPerElem][3]
>(stack.localResidual),
173 -stabilizationScaling );
174 stiffness.template upperBTDB< numNodesPerElem >( dNdX, -detJxW, stack.localJacobian );
177 template<
typename SUBREGION_TYPE,
178 typename CONSTITUTIVE_TYPE,
190 CONSTITUTIVE_TYPE::KernelWrapper::DiscretizationOps::template fillLowerBTDB< numNodesPerElem >( stack.localJacobian );
192 m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
193 for(
int localNode = 0; localNode < numSupportPoints; ++localNode )
195 for(
int dim = 0; dim < numDofPerTestSupportPoint; ++dim )
198 LvArray::integerConversion< localIndex >( stack.localRowDofIndex[ numDofPerTestSupportPoint * localNode + dim ] - m_dofRankOffset );
199 if( dof < 0 || dof >= m_matrix.numRows() )
201 m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( dof,
202 stack.localRowDofIndex,
203 stack.localJacobian[ numDofPerTestSupportPoint * localNode + dim ],
206 RAJA::atomicAdd< parallelDeviceAtomic >( &m_rhs[ dof ], stack.localResidual[ numDofPerTestSupportPoint * localNode + dim ] );
207 maxForce = fmax( maxForce, fabs( stack.localResidual[ numDofPerTestSupportPoint * localNode + dim ] ) );
213 template<
typename SUBREGION_TYPE,
214 typename CONSTITUTIVE_TYPE,
216 template<
typename POLICY,
217 typename KERNEL_TYPE >
220 KERNEL_TYPE
const & kernelComponent )
222 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.
double real64
64-bit floating point type.
GEOSX_GLOBALINDEX_TYPE globalIndex
Global index type (for indexing objects across MPI partitions).
GEOSX_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.