20 #ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_WELLS_SINGLEPHASEPERFORATIONFLUXLKERNELS_HPP
21 #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_WELLS_SINGLEPHASEPERFORATIONFLUXLKERNELS_HPP
23 #include "codingUtilities/Utilities.hpp"
25 #include "common/GEOS_RAJA_Interface.hpp"
28 #include "constitutive/fluid/singlefluid/SingleFluidLayouts.hpp"
30 #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp"
31 #include "constitutive/fluid/singlefluid/SingleFluidFields.hpp"
53 namespace isothermalSinglePhasePerforationFluxKernels
58 template<
integer IS_THERMAL >
72 fields::singlefluid::density,
73 fields::singlefluid::dDensity,
74 fields::singlefluid::viscosity,
75 fields::singlefluid::dViscosity >;
85 template<
typename VIEWTYPE >
90 constitutive::SingleFluidBase
const & fluid,
94 m_resPres( singlePhaseFlowAccessors.get( fields::flow::pressure {} )),
95 m_resDens( singleFluidAccessors.get( fields::singlefluid::density {} )),
96 m_dResDens( singleFluidAccessors.get( fields::singlefluid::dDensity {} )),
97 m_resVisc( singleFluidAccessors.get( fields::singlefluid::viscosity {} )),
98 m_dResVisc( singleFluidAccessors.get( fields::singlefluid::dViscosity {} )),
99 m_wellElemGravCoef( subRegion.getField< fields::well::gravityCoefficient >()),
100 m_wellElemPres( subRegion.getField< fields::well::pressure >()),
101 m_wellElemDens( fluid.density()),
102 m_dWellElemDens( fluid.dDensity()),
103 m_wellElemVisc( fluid.viscosity()),
104 m_dWellElemVisc( fluid.dViscosity()),
105 m_perfGravCoef( perforationData->getField< fields::well::gravityCoefficient >()),
106 m_perfWellElemIndex( perforationData->getField< fields::perforation::wellElementIndex >()),
107 m_perfTrans( perforationData->getField< fields::perforation::wellTransmissibility >()),
108 m_resElementRegion( perforationData->getField< fields::perforation::reservoirElementRegion >()),
109 m_resElementSubRegion( perforationData->getField< fields::perforation::reservoirElementSubRegion >()),
110 m_resElementIndex( perforationData->getField< fields::perforation::reservoirElementIndex >()),
111 m_perfRate( perforationData->getField< fields::well::perforationRate >()),
112 m_dPerfRate( perforationData->getField< fields::well::dPerforationRate >())
115 template<
typename FUNC = NoOpFunc >
119 computeFlux(
localIndex const iperf, FUNC && fluxKernelOp= NoOpFunc {} )
const
123 localIndex const er = m_resElementRegion[iperf];
124 localIndex const esr = m_resElementSubRegion[iperf];
125 localIndex const ei = m_resElementIndex[iperf];
128 localIndex const iwelem = m_perfWellElemIndex[iperf];
130 using Deriv = constitutive::singlefluid::DerivativeOffset;
131 using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< IS_THERMAL >;
135 real64 dPressure[2][DerivOffset::nDer]{};
139 m_perfRate[iperf] = 0.0;
140 for(
integer i = 0; i < 2; ++i )
142 for(
integer j=0; j<DerivOffset::nDer; j++ )
144 m_dPerfRate[iperf][i][j] = 0.0;
150 pressure[TAG::RES] = m_resPres[er][esr][ei];
151 dPressure[TAG::RES][DerivOffset::dP] = 1.0;
156 multiplier[TAG::RES] = 1;
161 pressure[TAG::WELL] = m_wellElemPres[iwelem];
162 dPressure[TAG::WELL][DerivOffset::dP] = 1.0;
164 real64 const gravD = ( m_perfGravCoef[iperf] - m_wellElemGravCoef[iwelem] );
165 pressure[TAG::WELL] += m_wellElemDens[iwelem][0] * gravD;
166 dPressure[TAG::WELL][DerivOffset::dP] += m_dWellElemDens[iwelem][0][Deriv::dP] * gravD;
167 if constexpr (IS_THERMAL)
169 dPressure[TAG::WELL][DerivOffset::dT] = m_dWellElemDens[iwelem][0][DerivOffset::dT] * gravD;
173 multiplier[TAG::WELL] = -1;
179 potDif += multiplier[i] * m_perfTrans[iperf] * pressure[i];
180 for(
integer j=0; j<DerivOffset::nDer; j++ )
182 m_dPerfRate[iperf][i][j] = multiplier[i] * m_perfTrans[iperf] * dPressure[i][j];
187 localIndex const k_up = (potDif >= 0) ? TAG::RES : TAG::WELL;
192 real64 dDensityUp[DerivOffset::nDer]{};
193 real64 dViscosityUp[DerivOffset::nDer]{};
194 for(
integer j=0; j<DerivOffset::nDer; j++ )
201 if( k_up == TAG::RES )
203 densityUp = m_resDens[er][esr][ei][0];
204 viscosityUp = m_resVisc[er][esr][ei][0];
205 for(
integer j=0; j<DerivOffset::nDer; j++ )
207 dDensityUp[j] = m_dResDens[er][esr][ei][0][j];
208 dViscosityUp[j] = m_dResVisc[er][esr][ei][0][j];
213 densityUp = m_wellElemDens[iwelem][0];
214 viscosityUp = m_wellElemVisc[iwelem][0];
215 for(
integer j=0; j<DerivOffset::nDer; j++ )
217 dDensityUp[j] = m_dWellElemDens[iwelem][0][j];
218 dViscosityUp[j] = m_dWellElemVisc[iwelem][0][j];
223 real64 const mobilityUp = densityUp / viscosityUp;
224 real64 dMobilityUp[DerivOffset::nDer]{};
225 for(
integer j=0; j<DerivOffset::nDer; j++ )
227 dMobilityUp[j] = dDensityUp[j] / viscosityUp
228 - mobilityUp / viscosityUp * dViscosityUp[j];
231 m_perfRate[iperf] = mobilityUp * potDif;
234 for(
integer j=0; j<DerivOffset::nDer; j++ )
236 m_dPerfRate[iperf][ke][j] *= mobilityUp;
239 for(
integer j=0; j<DerivOffset::nDer; j++ )
241 m_dPerfRate[iperf][k_up][j] += dMobilityUp[j]* potDif;
243 if constexpr ( IS_THERMAL )
245 fluxKernelOp( iwelem, er, esr, ei, potDif, m_perfRate[iperf], m_dPerfRate[iperf] );
255 template<
typename POLICY,
typename KERNEL_TYPE >
258 KERNEL_TYPE
const & kernelComponent )
263 kernelComponent.computeFlux( iperf );
269 ElementViewConst< arrayView1d< real64 const > >
const m_resPres;
270 ElementViewConst< arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > >
const m_resDens;
271 ElementViewConst< arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > >
const m_dResDens;
272 ElementViewConst< arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > >
const m_resVisc;
273 ElementViewConst< arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > >
const m_dResVisc;
309 template<
typename POLICY >
314 constitutive::SingleFluidBase
const & fluid,
317 integer constexpr IS_THERMAL = 0;
319 typename kernelType::SinglePhaseFlowAccessors singlePhaseFlowAccessors( elemManager, flowSolverName );
320 typename kernelType::SingleFluidAccessors singleFluidAccessors( elemManager, flowSolverName );
321 kernelType kernel( perforationData, subRegion, fluid, singlePhaseFlowAccessors, singleFluidAccessors );
322 kernelType::template launch< POLICY >( perforationData->
size(), kernel );
328 namespace thermalSinglePhasePerforationFluxKernels
331 using namespace constitutive;
335 template<
integer IS_THERMAL >
345 static constexpr
integer isThermal = IS_THERMAL;
354 fields::singlefluid::enthalpy,
355 fields::singlefluid::dEnthalpy >;
365 template<
typename VIEWTYPE >
370 SingleFluidBase
const & fluid,
375 :
Base( perforationData,
378 singlePhaseFlowAccessors,
379 singleFluidAccessors ),
380 m_wellElemEnthalpy( fluid.enthalpy()),
381 m_dWellElemEnthalpy( fluid.dEnthalpy()),
382 m_energyPerfFlux( perforationData->getField< fields::well::energyPerforationFlux >()),
383 m_dEnergyPerfFlux( perforationData->getField< fields::well::dEnergyPerforationFlux >()),
384 m_temp( thermalSinglePhaseFlowAccessors.get( fields::flow::temperature {} ) ),
385 m_resEnthalpy( thermalSingleFluidAccessors.get( fields::singlefluid::enthalpy {} ) ),
386 m_dResEnthalpy( thermalSingleFluidAccessors.get( fields::singlefluid::dEnthalpy {} ) )
390 template<
typename FUNC = NoOpFunc >
396 using Deriv = constitutive::singlefluid::DerivativeOffset;
397 using CP_Deriv =constitutive::singlefluid::DerivativeOffsetC< IS_THERMAL >;
399 m_energyPerfFlux[iperf]=0;
400 for(
integer ke = 0; ke < 2; ++ke )
402 for(
integer i = 0; i < CP_Deriv::nDer; ++i )
404 m_dEnergyPerfFlux[iperf][ke][i]=0;
409 real64 const potDiff,
real64 const flux, arraySlice2d< real64 >
const dFlux )
414 real64 const res_enthalpy = m_resEnthalpy[er][esr][ei][0];
416 m_energyPerfFlux[iperf] += flux * res_enthalpy;
419 m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dP] += dFlux[TAG::RES][CP_Deriv::dP] * res_enthalpy +
420 flux * m_dResEnthalpy[er][esr][ei][0][Deriv::dP];
421 m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dT] += dFlux[TAG::RES][CP_Deriv::dT] * res_enthalpy +
422 flux * m_dResEnthalpy[er][esr][ei][0][Deriv::dT];
424 m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dP] += dFlux[TAG::WELL][CP_Deriv::dP] * res_enthalpy;
425 m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dT] += dFlux[TAG::WELL][CP_Deriv::dT] * res_enthalpy;
430 real64 const wellelem_enthalpy = m_wellElemEnthalpy[iwelem][0];
431 m_energyPerfFlux[iperf] += flux * wellelem_enthalpy;
434 m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dP] += dFlux[TAG::RES][CP_Deriv::dP] * wellelem_enthalpy;
435 m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dT] += dFlux[TAG::RES][CP_Deriv::dT] * wellelem_enthalpy;
437 m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dP] += dFlux[TAG::WELL][CP_Deriv::dP] * wellelem_enthalpy
438 + flux * m_dWellElemEnthalpy[iwelem][0][Deriv::dP];
439 m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dT] += dFlux[TAG::WELL][CP_Deriv::dT] * wellelem_enthalpy
440 + flux * m_dWellElemEnthalpy[iwelem][0][Deriv::dT];
456 template<
typename POLICY,
typename KERNEL_TYPE >
459 KERNEL_TYPE
const & kernelComponent )
464 kernelComponent.computeFlux( iperf );
505 template<
typename POLICY >
510 SingleFluidBase
const & fluid,
513 integer constexpr IS_THERMAL = 1;
515 typename kernelType::SinglePhaseFlowAccessors singlePhaseFlowAccessors( elemManager, flowSolverName );
516 typename kernelType::SingleFluidAccessors singleFluidAccessors( elemManager, flowSolverName );
517 typename kernelType::ThermalSinglePhaseFlowAccessors thermalSinglePhaseFlowAccessors( elemManager, flowSolverName );
518 typename kernelType::ThermalSingleFluidAccessors thermalSingleFluidAccessors( elemManager, flowSolverName );
519 kernelType kernel( perforationData, subRegion, fluid, singlePhaseFlowAccessors, singleFluidAccessors, thermalSinglePhaseFlowAccessors, thermalSingleFluidAccessors );
520 kernelType::template launch< POLICY >( perforationData->
size(), kernel );
#define GEOS_HOST_DEVICE
Marks a host-device function.
#define GEOS_MARK_FUNCTION
Mark function with both Caliper and NVTX if enabled.
The ElementRegionManager class provides an interface to ObjectManagerBase in order to manage ElementR...
typename ElementViewAccessor< VIEWTYPE >::NestedViewTypeConst ElementViewConst
The ElementViewAccessor at the ElementRegionManager level is the type resulting from ElementViewAcces...
A struct to automatically construct and store element view accessors.
A struct to automatically construct and store element view accessors.
localIndex size() const
Get the "size" of the group, which determines the number of elements in resizable wrappers.
static void createAndLaunch(string const flowSolverName, PerforationData *const perforationData, ElementSubRegionBase const &subRegion, constitutive::SingleFluidBase const &fluid, ElementRegionManager &elemManager)
Create a new kernel and launch.
ElementRegionManager::ElementViewConst< VIEWTYPE > ElementViewConst
The type for element-based non-constitutive data parameters. Consists entirely of ArrayView's.
static constexpr integer isThermal
Compile time value for thermal option.
static void launch(localIndex const numElements, KERNEL_TYPE const &kernelComponent)
Performs the kernel launch.
static void createAndLaunch(string const flowSolverName, PerforationData *const perforationData, ElementSubRegionBase const &subRegion, SingleFluidBase const &fluid, ElementRegionManager &elemManager)
Create a new kernel and launch.
ElementViewConst< arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > > const m_resEnthalpy
Views on phase enthalpies.
static void launch(localIndex const numElements, KERNEL_TYPE const &kernelComponent)
Performs the kernel launch.
arrayView1d< real64 > const m_energyPerfFlux
Views on energy flux.
ElementRegionManager::ElementViewConst< VIEWTYPE > ElementViewConst
The type for element-based non-constitutive data parameters. Consists entirely of ArrayView's.
arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_wellElemEnthalpy
Views on well element properties.
ElementViewConst< arrayView1d< real64 const > > const m_temp
Views on temperature.
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
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.
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
ArrayView< T, 3, USD > arrayView3d
Alias for 3D array view.