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 >()),
113 m_perfStatus( perforationData->getField< fields::perforation::perforationStatus >())
116 template<
typename FUNC = NoOpFunc >
120 computeFlux(
localIndex const iperf, FUNC && fluxKernelOp= NoOpFunc {} )
const
123 using Deriv = constitutive::singlefluid::DerivativeOffset;
124 using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< IS_THERMAL >;
126 m_perfRate[iperf] = 0.0;
127 for(
integer i = 0; i < 2; ++i )
129 for(
integer j=0; j<DerivOffset::nDer; j++ )
131 m_dPerfRate[iperf][i][j] = 0.0;
135 if( !m_perfStatus[iperf] )
140 localIndex const er = m_resElementRegion[iperf];
141 localIndex const esr = m_resElementSubRegion[iperf];
142 localIndex const ei = m_resElementIndex[iperf];
145 localIndex const iwelem = m_perfWellElemIndex[iperf];
149 real64 dPressure[2][DerivOffset::nDer]{};
155 pressure[TAG::RES] = m_resPres[er][esr][ei];
156 dPressure[TAG::RES][DerivOffset::dP] = 1.0;
161 multiplier[TAG::RES] = 1;
166 pressure[TAG::WELL] = m_wellElemPres[iwelem];
167 dPressure[TAG::WELL][DerivOffset::dP] = 1.0;
169 real64 const gravD = ( m_perfGravCoef[iperf] - m_wellElemGravCoef[iwelem] );
170 pressure[TAG::WELL] += m_wellElemDens[iwelem][0] * gravD;
171 dPressure[TAG::WELL][DerivOffset::dP] += m_dWellElemDens[iwelem][0][Deriv::dP] * gravD;
172 if constexpr (IS_THERMAL)
174 dPressure[TAG::WELL][DerivOffset::dT] = m_dWellElemDens[iwelem][0][DerivOffset::dT] * gravD;
178 multiplier[TAG::WELL] = -1;
184 potDif += multiplier[i] * m_perfTrans[iperf] * pressure[i];
185 for(
integer j=0; j<DerivOffset::nDer; j++ )
187 m_dPerfRate[iperf][i][j] = multiplier[i] * m_perfTrans[iperf] * dPressure[i][j];
192 localIndex const k_up = (potDif >= 0) ? TAG::RES : TAG::WELL;
197 real64 dDensityUp[DerivOffset::nDer]{};
198 real64 dViscosityUp[DerivOffset::nDer]{};
199 for(
integer j=0; j<DerivOffset::nDer; j++ )
206 if( k_up == TAG::RES )
208 densityUp = m_resDens[er][esr][ei][0];
209 viscosityUp = m_resVisc[er][esr][ei][0];
210 for(
integer j=0; j<DerivOffset::nDer; j++ )
212 dDensityUp[j] = m_dResDens[er][esr][ei][0][j];
213 dViscosityUp[j] = m_dResVisc[er][esr][ei][0][j];
218 densityUp = m_wellElemDens[iwelem][0];
219 viscosityUp = m_wellElemVisc[iwelem][0];
220 for(
integer j=0; j<DerivOffset::nDer; j++ )
222 dDensityUp[j] = m_dWellElemDens[iwelem][0][j];
223 dViscosityUp[j] = m_dWellElemVisc[iwelem][0][j];
228 real64 const mobilityUp = densityUp / viscosityUp;
229 real64 dMobilityUp[DerivOffset::nDer]{};
230 for(
integer j=0; j<DerivOffset::nDer; j++ )
232 dMobilityUp[j] = dDensityUp[j] / viscosityUp
233 - mobilityUp / viscosityUp * dViscosityUp[j];
236 m_perfRate[iperf] = mobilityUp * potDif;
239 for(
integer j=0; j<DerivOffset::nDer; j++ )
241 m_dPerfRate[iperf][ke][j] *= mobilityUp;
244 for(
integer j=0; j<DerivOffset::nDer; j++ )
246 m_dPerfRate[iperf][k_up][j] += dMobilityUp[j]* potDif;
248 if constexpr ( IS_THERMAL )
250 fluxKernelOp( iwelem, er, esr, ei, potDif, m_perfRate[iperf], m_dPerfRate[iperf] );
260 template<
typename POLICY,
typename KERNEL_TYPE >
263 KERNEL_TYPE
const & kernelComponent )
268 kernelComponent.computeFlux( iperf );
274 ElementViewConst< arrayView1d< real64 const > >
const m_resPres;
275 ElementViewConst< arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > >
const m_resDens;
276 ElementViewConst< arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > >
const m_dResDens;
277 ElementViewConst< arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > >
const m_resVisc;
278 ElementViewConst< arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > >
const m_dResVisc;
315 template<
typename POLICY >
320 constitutive::SingleFluidBase
const & fluid,
323 integer constexpr IS_THERMAL = 0;
325 typename kernelType::SinglePhaseFlowAccessors singlePhaseFlowAccessors( elemManager, flowSolverName );
326 typename kernelType::SingleFluidAccessors singleFluidAccessors( elemManager, flowSolverName );
327 kernelType kernel( perforationData, subRegion, fluid, singlePhaseFlowAccessors, singleFluidAccessors );
328 kernelType::template launch< POLICY >( perforationData->
size(), kernel );
334 namespace thermalSinglePhasePerforationFluxKernels
337 using namespace constitutive;
341 template<
integer IS_THERMAL >
351 static constexpr
integer isThermal = IS_THERMAL;
360 fields::singlefluid::enthalpy,
361 fields::singlefluid::dEnthalpy >;
371 template<
typename VIEWTYPE >
376 SingleFluidBase
const & fluid,
381 :
Base( perforationData,
384 singlePhaseFlowAccessors,
385 singleFluidAccessors ),
386 m_wellElemEnthalpy( fluid.enthalpy()),
387 m_dWellElemEnthalpy( fluid.dEnthalpy()),
388 m_energyPerfFlux( perforationData->getField< fields::well::energyPerforationFlux >()),
389 m_dEnergyPerfFlux( perforationData->getField< fields::well::dEnergyPerforationFlux >()),
390 m_temp( thermalSinglePhaseFlowAccessors.get( fields::flow::temperature {} ) ),
391 m_resEnthalpy( thermalSingleFluidAccessors.get( fields::singlefluid::enthalpy {} ) ),
392 m_dResEnthalpy( thermalSingleFluidAccessors.get( fields::singlefluid::dEnthalpy {} ) )
396 template<
typename FUNC = NoOpFunc >
402 using Deriv = constitutive::singlefluid::DerivativeOffset;
403 using CP_Deriv =constitutive::singlefluid::DerivativeOffsetC< IS_THERMAL >;
405 m_energyPerfFlux[iperf]=0;
406 for(
integer ke = 0; ke < 2; ++ke )
408 for(
integer i = 0; i < CP_Deriv::nDer; ++i )
410 m_dEnergyPerfFlux[iperf][ke][i]=0;
415 real64 const potDiff,
real64 const flux, arraySlice2d< real64 >
const dFlux )
420 real64 const res_enthalpy = m_resEnthalpy[er][esr][ei][0];
422 m_energyPerfFlux[iperf] += flux * res_enthalpy;
425 m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dP] += dFlux[TAG::RES][CP_Deriv::dP] * res_enthalpy +
426 flux * m_dResEnthalpy[er][esr][ei][0][Deriv::dP];
427 m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dT] += dFlux[TAG::RES][CP_Deriv::dT] * res_enthalpy +
428 flux * m_dResEnthalpy[er][esr][ei][0][Deriv::dT];
430 m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dP] += dFlux[TAG::WELL][CP_Deriv::dP] * res_enthalpy;
431 m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dT] += dFlux[TAG::WELL][CP_Deriv::dT] * res_enthalpy;
436 real64 const wellelem_enthalpy = m_wellElemEnthalpy[iwelem][0];
437 m_energyPerfFlux[iperf] += flux * wellelem_enthalpy;
440 m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dP] += dFlux[TAG::RES][CP_Deriv::dP] * wellelem_enthalpy;
441 m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dT] += dFlux[TAG::RES][CP_Deriv::dT] * wellelem_enthalpy;
443 m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dP] += dFlux[TAG::WELL][CP_Deriv::dP] * wellelem_enthalpy
444 + flux * m_dWellElemEnthalpy[iwelem][0][Deriv::dP];
445 m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dT] += dFlux[TAG::WELL][CP_Deriv::dT] * wellelem_enthalpy
446 + flux * m_dWellElemEnthalpy[iwelem][0][Deriv::dT];
462 template<
typename POLICY,
typename KERNEL_TYPE >
465 KERNEL_TYPE
const & kernelComponent )
470 kernelComponent.computeFlux( iperf );
511 template<
typename POLICY >
516 SingleFluidBase
const & fluid,
519 integer constexpr IS_THERMAL = 1;
521 typename kernelType::SinglePhaseFlowAccessors singlePhaseFlowAccessors( elemManager, flowSolverName );
522 typename kernelType::SingleFluidAccessors singleFluidAccessors( elemManager, flowSolverName );
523 typename kernelType::ThermalSinglePhaseFlowAccessors thermalSinglePhaseFlowAccessors( elemManager, flowSolverName );
524 typename kernelType::ThermalSingleFluidAccessors thermalSingleFluidAccessors( elemManager, flowSolverName );
525 kernelType kernel( perforationData, subRegion, fluid, singlePhaseFlowAccessors, singleFluidAccessors, thermalSinglePhaseFlowAccessors, thermalSingleFluidAccessors );
526 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).
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
int integer
Signed integer type.
ArrayView< T, 3, USD > arrayView3d
Alias for 3D array view.