20 #ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_COUPLEDRESERVOIRANDWELLS_HPP
21 #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_COUPLEDRESERVOIRANDWELLS_HPP
25 #include "common/GEOS_RAJA_Interface.hpp"
26 #include "constitutive/fluid/multifluid/Layouts.hpp"
34 namespace coupledReservoirAndWellKernels
37 using namespace constitutive;
44 template<
integer NC,
integer IS_THERMAL >
51 static constexpr
integer resNumDOF = NC+1+IS_THERMAL;
60 using CP_Deriv = multifluid::DerivativeOffsetC< NC, IS_THERMAL >;
67 static constexpr
integer numDof = WJ_COFFSET::nDer;
70 static constexpr
integer numEqn = WJ_ROFFSET::nEqn - 2;
88 string const wellDofKey,
92 MultiFluidBase
const & fluid,
96 bool const & detectCrossflow,
97 integer & numCrossFlowPerforations,
98 BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags )
101 m_numPhases ( fluid.numFluidPhases()),
102 m_rankOffset( rankOffset ),
103 m_compPerfRate( perforationData->getField< fields::well::compPerforationRate >() ),
104 m_dCompPerfRate( perforationData->getField< fields::well::dCompPerforationRate >() ),
105 m_perfWellElemIndex( perforationData->getField< fields::perforation::wellElementIndex >() ),
106 m_perfStatus( perforationData->getField< fields::perforation::perforationStatus >() ),
107 m_wellElemDofNumber( subRegion.getReference<
array1d<
globalIndex > >( wellDofKey ) ),
108 m_resElemDofNumber( resDofNumber ),
109 m_resElementRegion( perforationData->getField< fields::perforation::reservoirElementRegion >() ),
110 m_resElementSubRegion( perforationData->getField< fields::perforation::reservoirElementSubRegion >() ),
111 m_resElementIndex( perforationData->getField< fields::perforation::reservoirElementIndex >() ),
112 m_localRhs( localRhs ),
113 m_localMatrix( localMatrix ),
114 m_detectCrossflow( detectCrossflow ),
115 m_numCrossFlowPerforations( numCrossFlowPerforations ),
116 m_useTotalMassEquation ( kernelFlags.isSet( isothermalCompositionalMultiphaseBaseKernels::KernelFlags::TotalMassEquation ) )
128 template<
typename FUNC = NoOpFunc >
132 FUNC && compFluxKernelOp = NoOpFunc{} )
const
135 using namespace compositionalMultiphaseUtilities;
136 if( m_perfStatus[iperf ] )
146 localIndex const er = m_resElementRegion[iperf];
147 localIndex const esr = m_resElementSubRegion[iperf];
148 localIndex const ei = m_resElementIndex[iperf];
151 localIndex const iwelem = m_perfWellElemIndex[iperf];
152 globalIndex const resOffset = m_resElemDofNumber[er][esr][ei];
153 globalIndex const wellElemOffset = m_wellElemDofNumber[iwelem];
155 for(
integer ic = 0; ic < numComp; ++ic )
157 eqnRowIndices[TAG::RES * numComp + ic] = LvArray::integerConversion< localIndex >( resOffset - m_rankOffset ) + ic;
158 eqnRowIndices[TAG::WELL * numComp + ic] = LvArray::integerConversion< localIndex >( wellElemOffset - m_rankOffset ) + WJ_ROFFSET::MASSBAL + ic;
161 for(
integer jdof = 0; jdof < NC+1; ++jdof )
163 dofColIndices[TAG::RES * resNumDOF + jdof] = resOffset + jdof;
164 dofColIndices[TAG::WELL * resNumDOF + jdof] = wellElemOffset + WJ_COFFSET::dP + jdof;
167 if constexpr ( IS_THERMAL )
169 dofColIndices[TAG::RES * resNumDOF + NC+1 ] = resOffset + NC+1;
170 dofColIndices[TAG::WELL * resNumDOF + NC+1 ] = wellElemOffset + WJ_COFFSET::dT;
173 for(
integer ic = 0; ic < numComp; ++ic )
175 localPerf[TAG::RES * numComp + ic] = m_dt * m_compPerfRate[iperf][ic];
176 localPerf[TAG::WELL * numComp + ic] = -m_dt * m_compPerfRate[iperf][ic];
178 if( m_detectCrossflow )
180 if( m_compPerfRate[iperf][ic] > LvArray::NumericLimits< real64 >::epsilon )
182 m_numCrossFlowPerforations += 1;
185 for(
integer ke = 0; ke < 2; ++ke )
187 localIndex localDofIndexPres = ke * resNumDOF;
189 localPerfJacobian[TAG::RES * numComp + ic][localDofIndexPres] = m_dt * m_dCompPerfRate[iperf][ke][ic][CP_Deriv::dP];
190 localPerfJacobian[TAG::WELL * numComp + ic][localDofIndexPres] = -m_dt * m_dCompPerfRate[iperf][ke][ic][CP_Deriv::dP];
191 for(
integer jc = 0; jc < numComp; ++jc )
193 localIndex const localDofIndexComp = localDofIndexPres + jc + 1;
195 localPerfJacobian[TAG::RES * numComp + ic][localDofIndexComp] = m_dt * m_dCompPerfRate[iperf][ke][ic][CP_Deriv::dC+jc];
196 localPerfJacobian[TAG::WELL * numComp + ic][localDofIndexComp] = -m_dt * m_dCompPerfRate[iperf][ke][ic][CP_Deriv::dC+jc];
198 if constexpr ( IS_THERMAL )
200 localIndex localDofIndexTemp = localDofIndexPres + NC + 1;
201 localPerfJacobian[TAG::RES * numComp + ic][localDofIndexTemp] = m_dt * m_dCompPerfRate[iperf][ke][ic][CP_Deriv::dT];
202 localPerfJacobian[TAG::WELL * numComp + ic][localDofIndexTemp] = -m_dt * m_dCompPerfRate[iperf][ke][ic][CP_Deriv::dT];
207 if( m_useTotalMassEquation )
210 stackArray1d< real64, 2 * resNumDOF > work( 2 * resNumDOF );
215 for(
localIndex i = 0; i < localPerf.size(); ++i )
217 if( eqnRowIndices[i] >= 0 && eqnRowIndices[i] < m_localMatrix.numRows() )
219 m_localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >( eqnRowIndices[i],
220 dofColIndices.data(),
221 localPerfJacobian[i].dataIfContiguous(),
223 RAJA::atomicAdd( parallelDeviceAtomic{}, &m_localRhs[eqnRowIndices[i]], localPerf[i] );
226 compFluxKernelOp( resOffset, wellElemOffset, dofColIndices, iwelem );
239 template<
typename POLICY,
typename KERNEL_TYPE >
242 KERNEL_TYPE
const & kernelComponent )
247 kernelComponent.computeFlux( ie );
277 bool const m_detectCrossflow;
278 integer & m_numCrossFlowPerforations;
279 integer const m_useTotalMassEquation;
302 template<
typename POLICY >
307 string const wellDofKey,
311 MultiFluidBase
const & fluid,
312 BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags,
313 bool const & detectCrossflow,
314 integer & numCrossFlowPerforations,
319 isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComps, [&](
auto NC )
321 integer constexpr NUM_COMP = NC();
324 kernelType kernel( dt, rankOffset, wellDofKey, subRegion, resDofNumber, perforationData,
325 fluid, localRhs, localMatrix, detectCrossflow, numCrossFlowPerforations, kernelFlags );
326 kernelType::template launch< POLICY >( perforationData->
size(), kernel );
338 template<
integer NC,
integer IS_THERMAL >
345 static constexpr
integer resNumDOF = NC+1+IS_THERMAL;
354 using CP_Deriv = multifluid::DerivativeOffsetC< NC, IS_THERMAL >;
359 using Base::m_localRhs;
360 using Base::m_localMatrix;
361 using Base::m_rankOffset;
366 static constexpr
integer numDof = WJ_COFFSET::nDer;
369 static constexpr
integer numEqn = WJ_ROFFSET::nEqn - 2;
388 string const wellDofKey,
392 MultiFluidBase
const & fluid,
395 bool const & detectCrossflow,
396 integer & numCrossFlowPerforations,
397 BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags )
408 numCrossFlowPerforations,
410 m_isProducer( isProducer ),
411 m_globalWellElementIndex( subRegion.getGlobalWellElementIndex() ),
412 m_energyPerfFlux( perforationData->getField< fields::well::energyPerforationFlux >()),
413 m_dEnergyPerfFlux( perforationData->getField< fields::well::dEnergyPerforationFlux >())
430 Base::computeFlux( iperf, [&] (
globalIndex const & resOffset,
440 if( m_globalWellElementIndex[iwelem] == 0 )
451 eqnRowIndices[TAG::RES ] = LvArray::integerConversion< localIndex >( resOffset - m_rankOffset ) + NC + 1;
452 eqnRowIndices[TAG::WELL ] = LvArray::integerConversion< localIndex >( wellElemOffset - m_rankOffset ) + WJ_ROFFSET::ENERGYBAL;
455 localPerf[TAG::RES ] = m_dt * m_energyPerfFlux[iperf];
456 localPerf[TAG::WELL ] = -m_dt * m_energyPerfFlux[iperf];
458 for(
integer ke = 0; ke < 2; ++ke )
460 localIndex localDofIndexPres = ke * resNumDOF;
461 localPerfJacobian[TAG::RES ][localDofIndexPres] = m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dP];
462 localPerfJacobian[TAG::WELL ][localDofIndexPres] = -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dP];
465 for(
integer ic = 0; ic < numComp; ++ic )
467 localIndex const localDofIndexComp = localDofIndexPres + ic + 1;
468 localPerfJacobian[TAG::RES ][localDofIndexComp] = m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dC+ic];
469 localPerfJacobian[TAG::WELL][localDofIndexComp] = -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dC+ic];
471 localPerfJacobian[TAG::RES ][localDofIndexPres+NC+1] = m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dT];
472 localPerfJacobian[TAG::WELL][localDofIndexPres+NC+1] = -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dT];
476 for(
localIndex i = 0; i < localPerf.size(); ++i )
478 if( eqnRowIndices[i] >= 0 && eqnRowIndices[i] < m_localMatrix.numRows() )
480 m_localMatrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( eqnRowIndices[i],
481 dofColIndices.data(),
482 localPerfJacobian[i].dataIfContiguous(),
484 RAJA::atomicAdd( parallelDeviceAtomic{}, &m_localRhs[eqnRowIndices[i]], localPerf[i] );
499 template<
typename POLICY,
typename KERNEL_TYPE >
502 KERNEL_TYPE
const & kernelComponent )
507 kernelComponent.computeFlux( ie );
545 template<
typename POLICY >
551 string const wellDofKey,
555 MultiFluidBase
const & fluid,
556 BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags,
557 bool const & detectCrossflow,
558 integer & numCrossFlowPerforations,
563 isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComps, [&](
auto NC )
565 integer constexpr NUM_COMP = NC();
568 kernelType kernel( dt, isProducer, rankOffset, wellDofKey, subRegion, resDofNumber, perforationData,
569 fluid, localRhs, localMatrix, detectCrossflow, numCrossFlowPerforations, kernelFlags );
570 kernelType::template launch< POLICY >( perforationData->
size(), kernel );
GEOS_HOST_DEVICE void shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum(integer const numRowsToShift, integer const numRowsInBlock, integer const numColsInBlock, integer const numBlocks, MATRIX &&mat, VEC &&work)
In each block, shift the elements from 0 to numRowsToShift-1, shifts all rows one position ahead and ...
GEOS_HOST_DEVICE void shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum(integer const numRowsToShift, integer const numRowsInBlock, integer const numBlocks, VEC &&v)
In each block, shift the elements from 0 to numRowsToShift-1 one position ahead and replaces the firs...
#define GEOS_HOST_DEVICE
Marks a host-device function.
#define GEOS_MARK_FUNCTION
Mark function with both Caliper and NVTX if enabled.
typename ElementViewAccessor< VIEWTYPE >::NestedViewTypeConst ElementViewConst
The ElementViewAccessor at the ElementRegionManager level is the type resulting from ElementViewAcces...
This class describes a collection of local well elements and perforations.
static void createAndLaunch(integer const numComps, real64 const dt, globalIndex const rankOffset, string const wellDofKey, WellElementSubRegion const &subRegion, ElementRegionManager::ElementViewConst< arrayView1d< globalIndex const > > const resDofNumber, PerforationData const *const perforationData, MultiFluidBase const &fluid, BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags, bool const &detectCrossflow, integer &numCrossFlowPerforations, arrayView1d< real64 > const &localRhs, CRSMatrixView< real64, globalIndex const > const &localMatrix)
Create a new kernel and launch.
integer const m_numPhases
Number of phases.
GEOS_HOST_DEVICE void computeFlux(localIndex const iperf, FUNC &&compFluxKernelOp=NoOpFunc{}) const
Compute the local flux contributions to the residual and Jacobian.
IsothermalCompositionalMultiPhaseFluxKernel(real64 const dt, globalIndex const rankOffset, string const wellDofKey, WellElementSubRegion const &subRegion, ElementRegionManager::ElementViewConst< arrayView1d< globalIndex const > > const resDofNumber, PerforationData const *const perforationData, MultiFluidBase const &fluid, arrayView1d< real64 > const &localRhs, CRSMatrixView< real64, globalIndex const > const &localMatrix, bool const &detectCrossflow, integer &numCrossFlowPerforations, BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags)
Constructor for the kernel interface.
real64 const m_dt
Time step size.
static void launch(localIndex const numElements, KERNEL_TYPE const &kernelComponent)
Performs the kernel launch.
static void createAndLaunch(integer const numComps, integer const isProducer, real64 const dt, globalIndex const rankOffset, string const wellDofKey, WellElementSubRegion const &subRegion, ElementRegionManager::ElementViewConst< arrayView1d< globalIndex const > > const resDofNumber, PerforationData const *const perforationData, MultiFluidBase const &fluid, BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags, bool const &detectCrossflow, integer &numCrossFlowPerforations, arrayView1d< real64 > const &localRhs, CRSMatrixView< real64, globalIndex const > const &localMatrix)
Create a new kernel and launch.
GEOS_HOST_DEVICE void computeFlux(localIndex const iperf) const
Compute the local flux contributions to the residual and Jacobian.
ThermalCompositionalMultiPhaseFluxKernel(real64 const dt, integer const isProducer, globalIndex const rankOffset, string const wellDofKey, WellElementSubRegion const &subRegion, ElementRegionManager::ElementViewConst< arrayView1d< globalIndex const > > const resDofNumber, PerforationData const *const perforationData, MultiFluidBase const &fluid, arrayView1d< real64 > const &localRhs, CRSMatrixView< real64, globalIndex const > const &localMatrix, bool const &detectCrossflow, integer &numCrossFlowPerforations, BitFlags< isothermalCompositionalMultiphaseBaseKernels::KernelFlags > kernelFlags)
Constructor for the kernel interface.
arrayView1d< real64 const > const m_energyPerfFlux
Views on energy flux.
static void launch(localIndex const numElements, KERNEL_TYPE const &kernelComponent)
Performs the kernel launch.
integer const m_isProducer
Well type.
arrayView1d< globalIndex const > m_globalWellElementIndex
Global index of local element.
localIndex size() const
Get the "size" of the group, which determines the number of elements in resizable wrappers.
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
StackArray< T, 2, MAXSIZE > stackArray2d
Alias for 2D stack array.
GEOS_GLOBALINDEX_TYPE globalIndex
Global index type (for indexing objects across MPI partitions).
StackArray< T, 1, MAXSIZE > stackArray1d
Alias for 1D stack array.
double real64
64-bit floating point type.
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
LvArray::CRSMatrixView< T, COL_INDEX, INDEX_TYPE const, LvArray::ChaiBuffer > CRSMatrixView
Alias for CRS Matrix View.
ArrayView< T, 4, USD > arrayView4d
Alias for 4D array view.
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
int integer
Signed integer type.
Array< T, 1 > array1d
Alias for 1D array.
ArrayView< T, 3, USD > arrayView3d
Alias for 3D array view.