20 #ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_ACCUMULATIONKERNEL_HPP
21 #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONAL_ACCUMULATIONKERNEL_HPP
23 #include "codingUtilities/Utilities.hpp"
26 #include "common/GEOS_RAJA_Interface.hpp"
27 #include "constitutive/solid/CoupledSolidBase.hpp"
28 #include "constitutive/fluid/multifluid/MultiFluidBase.hpp"
33 #include "physicsSolvers/fluidFlow/kernels/compositional/KernelLaunchSelectors.hpp"
38 namespace isothermalCompositionalMultiphaseBaseKernels
41 static constexpr
real64 minDensForDivision = 1e-10;
43 enum class KernelFlags
45 SimpleAccumulation = 1 << 0,
46 TotalMassEquation = 1 << 1,
64 template<
integer NUM_COMP,
integer NUM_DOF >
93 constitutive::MultiFluidBase
const & fluid,
94 constitutive::CoupledSolidBase
const & solid,
97 BitFlags< KernelFlags >
const kernelFlags )
102 m_volume( subRegion.getElementVolume() ),
104 m_porosity_n( solid.getPorosity_n() ),
105 m_dPoro_dPres( solid.getDporosity_dPressure() ),
106 m_dCompFrac_dCompDens( subRegion.getField< fields::flow::dGlobalCompFraction_dGlobalCompDensity >() ),
107 m_phaseVolFrac( subRegion.getField< fields::flow::phaseVolumeFraction >() ),
108 m_dPhaseVolFrac( subRegion.getField< fields::flow::dPhaseVolumeFraction >() ),
110 m_dPhaseDens( fluid.dPhaseDensity() ),
112 m_dPhaseCompFrac( fluid.dPhaseCompFraction() ),
113 m_compDens( subRegion.getField< fields::flow::globalCompDensity >() ),
114 m_compAmount_n( subRegion.getField< fields::flow::compAmount_n >() ),
117 m_kernelFlags( kernelFlags )
190 template<
typename FUNC = NoOpFunc >
194 FUNC && phaseAmountKernelOp = NoOpFunc{} )
const
196 if( m_kernelFlags.isSet( KernelFlags::SimpleAccumulation ) )
203 real64 const compAmount_n = m_compAmount_n[ei][ic];
218 using Deriv = constitutive::multifluid::DerivativeOffset;
224 arraySlice2d<
real64 const, compflow::USD_PHASE_DC - 1 > dPhaseVolFrac = m_dPhaseVolFrac[ei];
227 arraySlice2d<
real64 const, constitutive::multifluid::USD_PHASE_DC - 2 > dPhaseDens = m_dPhaseDens[ei][0];
230 arraySlice3d<
real64 const, constitutive::multifluid::USD_PHASE_COMP_DC - 2 > dPhaseCompFrac = m_dPhaseCompFrac[ei][0];
245 real64 const phaseAmount = stack.
poreVolume * phaseVolFrac[ip] * phaseDens[ip];
248 + stack.
poreVolume * ( dPhaseVolFrac[ip][Deriv::dP] * phaseDens[ip]
249 + phaseVolFrac[ip] * dPhaseDens[ip][Deriv::dP] );
252 applyChainRule(
numComp, dCompFrac_dCompDens, dPhaseDens[ip], dPhaseAmount_dC, Deriv::dC );
255 dPhaseAmount_dC[jc] = dPhaseAmount_dC[jc] * phaseVolFrac[ip]
256 + phaseDens[ip] * dPhaseVolFrac[ip][Deriv::dC + jc];
264 real64 const phaseCompAmount = phaseAmount * phaseCompFrac[ip][ic];
266 real64 const dPhaseCompAmount_dP = dPhaseAmount_dP * phaseCompFrac[ip][ic]
267 + phaseAmount * dPhaseCompFrac[ip][ic][Deriv::dP];
276 applyChainRule(
numComp, dCompFrac_dCompDens, dPhaseCompFrac[ip][ic], dPhaseCompFrac_dC, Deriv::dC );
279 real64 const dPhaseCompAmount_dC = dPhaseCompFrac_dC[jc] * phaseAmount
280 + phaseCompFrac[ip][ic] * dPhaseAmount_dC[jc];
288 phaseAmountKernelOp( ip, phaseAmount, dPhaseAmount_dP, dPhaseAmount_dC );
301 template<
typename FUNC = NoOpFunc >
305 FUNC && phaseVolFractionSumKernelOp = NoOpFunc{} )
const
307 using Deriv = constitutive::multifluid::DerivativeOffset;
310 arraySlice2d<
real64 const, compflow::USD_PHASE_DC - 1 > dPhaseVolFrac = m_dPhaseVolFrac[ei];
312 real64 oneMinusPhaseVolFracSum = 1.0;
317 oneMinusPhaseVolFracSum -= phaseVolFrac[ip];
328 phaseVolFractionSumKernelOp( oneMinusPhaseVolFracSum );
348 using namespace compositionalMultiphaseUtilities;
350 if( m_kernelFlags.isSet( KernelFlags::TotalMassEquation ) )
363 for(
integer i = 0; i < numRows; ++i )
380 template<
typename POLICY,
typename KERNEL_TYPE >
383 KERNEL_TYPE
const & kernelComponent )
389 if( kernelComponent.elemGhostRank( ei ) >= 0 )
394 typename KERNEL_TYPE::StackVariables stack;
396 kernelComponent.setup( ei, stack );
397 kernelComponent.computeAccumulation( ei, stack );
398 kernelComponent.computeVolumeBalance( ei, stack );
399 kernelComponent.complete( ei, stack );
451 BitFlags< KernelFlags >
const m_kernelFlags;
474 template<
typename POLICY >
479 BitFlags< KernelFlags > kernelFlags,
482 constitutive::MultiFluidBase
const & fluid,
483 constitutive::CoupledSolidBase
const & solid,
487 internal::kernelLaunchSelectorCompSwitch( numComps, [&] (
auto NC )
489 integer constexpr NUM_COMP = NC();
490 integer constexpr NUM_DOF = NC()+1;
493 fluid, solid, localMatrix, localRhs, kernelFlags );
#define GEOS_HOST_DEVICE
Marks a host-device function.
#define GEOS_UNUSED_PARAM(X)
Mark an unused argument and silence compiler warnings.
#define GEOS_MARK_FUNCTION
Mark function with both Caliper and NVTX if enabled.
localIndex size() const
Get the "size" of the group, which determines the number of elements in resizable wrappers.
static void createAndLaunch(integer const numComps, integer const numPhases, globalIndex const rankOffset, BitFlags< KernelFlags > kernelFlags, string const dofKey, ElementSubRegionBase const &subRegion, constitutive::MultiFluidBase const &fluid, constitutive::CoupledSolidBase const &solid, CRSMatrixView< real64, globalIndex const > const &localMatrix, arrayView1d< real64 > const &localRhs)
Create a new kernel and launch.
Define the interface for the assembly kernel in charge of accumulation and volume balance.
static constexpr integer numDof
Compute time value for the number of degrees of freedom.
arrayView2d< real64 const, compflow::USD_PHASE > const m_phaseVolFrac
Views on the phase volume fractions.
GEOS_HOST_DEVICE void complete(localIndex const GEOS_UNUSED_PARAM(ei), StackVariables &stack) const
Performs the complete phase for the kernel.
arrayView1d< integer const > const m_elemGhostRank
View on the ghost ranks.
CRSMatrixView< real64, globalIndex const > const m_localMatrix
View on the local CRS matrix.
arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > const m_phaseCompFrac
Views on the phase component fraction.
static void launch(localIndex const numElems, KERNEL_TYPE const &kernelComponent)
Performs the kernel launch.
static constexpr integer numEqn
Compute time value for the number of equations.
AccumulationKernel(localIndex const numPhases, globalIndex const rankOffset, string const dofKey, ElementSubRegionBase const &subRegion, constitutive::MultiFluidBase const &fluid, constitutive::CoupledSolidBase const &solid, CRSMatrixView< real64, globalIndex const > const &localMatrix, arrayView1d< real64 > const &localRhs, BitFlags< KernelFlags > const kernelFlags)
Constructor.
GEOS_HOST_DEVICE void computeVolumeBalance(localIndex const ei, StackVariables &stack, FUNC &&phaseVolFractionSumKernelOp=NoOpFunc{}) const
Compute the local volume balance contributions to the residual and Jacobian.
arrayView3d< real64 const, compflow::USD_COMP_DC > const m_dCompFrac_dCompDens
Views on the derivatives of comp fractions wrt component density.
arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > const m_phaseDens
Views on the phase densities.
arrayView2d< real64 const > const m_porosity
Views on the porosity.
arrayView1d< real64 > const m_localRhs
View on the local RHS.
GEOS_HOST_DEVICE void computeAccumulation(localIndex const ei, StackVariables &stack, FUNC &&phaseAmountKernelOp=NoOpFunc{}) const
Compute the local accumulation contributions to the residual and Jacobian.
GEOS_HOST_DEVICE integer elemGhostRank(localIndex const ei) const
Getter for the ghost rank of an element.
GEOS_HOST_DEVICE void setup(localIndex const ei, StackVariables &stack) const
Performs the setup phase for the kernel.
arrayView1d< real64 const > const m_volume
View on the element volumes.
integer const m_numPhases
Number of fluid phases.
static constexpr integer numComp
Compile time value for the number of components.
globalIndex const m_rankOffset
Offset for my MPI rank.
arrayView1d< globalIndex const > const m_dofNumber
View on the dof numbers.
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).
ArraySlice< T, 2, USD > arraySlice2d
Alias for 2D array slice.
ArraySlice< T, 3, USD > arraySlice3d
Alias for 3D array slice.
ArrayView< T, 5, USD > arrayView5d
Alias for 5D array view.
double real64
64-bit floating point type.
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
ArraySlice< T, 1, USD > arraySlice1d
Alias for 1D array slice.
std::int32_t integer
Signed integer type.
ArrayView< T, 4, USD > arrayView4d
Alias for 4D array view.
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
Array< T, 1 > array1d
Alias for 1D array.
ArrayView< T, 3, USD > arrayView3d
Alias for 3D array view.
Kernel variables (dof numbers, jacobian and residual) located on the stack.
globalIndex dofIndices[numDof]
Indices of the matrix rows/columns corresponding to the dofs in this element.
real64 dPoreVolume_dPres
Derivative of pore volume with respect to pressure.
localIndex localRow
Index of the local row corresponding to this element.
real64 poreVolume
Pore volume at time n+1.
real64 localJacobian[numEqn][numDof]
C-array storage for the element local Jacobian matrix (all equations except volume balance,...
real64 localResidual[numEqn]
C-array storage for the element local residual vector (all equations except volume balance)