Skip to content

feat: enable crossflow logic for producers #3800

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Sep 9, 2025
Merged

feat: enable crossflow logic for producers #3800

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
7a73be1
feat: enable crossflow logic for producers
Aug 29, 2025
01e4070
Merge branch 'develop' into pt/prod-crossflow
paveltomin Sep 9, 2025
39eb083
Merge branch 'develop' into pt/prod-crossflow
paveltomin Sep 9, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 4 additions & 5 deletions src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1709,9 +1709,6 @@ void CompositionalMultiphaseWell::computePerforationRates( real64 const & time_n
WellControls const & wellControls = getWellControls( subRegion );
if( wellControls.getWellStatus() == WellControls::Status::OPEN && !m_keepVariablesConstantDuringInitStep )
{

bool const disableReservoirToWellFlow = wellControls.isInjector() and !wellControls.isCrossflowEnabled();

string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() );
MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName );
bool const isThermal = fluid.isThermal();
Expand All @@ -1727,7 +1724,8 @@ void CompositionalMultiphaseWell::computePerforationRates( real64 const & time_n
subRegion,
fluid,
elemManager,
disableReservoirToWellFlow );
wellControls.isInjector(),
wellControls.isCrossflowEnabled() );
}
else
{
Expand All @@ -1739,7 +1737,8 @@ void CompositionalMultiphaseWell::computePerforationRates( real64 const & time_n
perforationData,
subRegion,
elemManager,
disableReservoirToWellFlow );
wellControls.isInjector(),
wellControls.isCrossflowEnabled() );
}
}
else
Expand Down
6 changes: 0 additions & 6 deletions src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -376,12 +376,6 @@ void WellControls::postInputInitialization()
<< EnumStrings< Control >::toString( Control::MASSRATE ),
InputError );

// 7) Make sure that the flag disabling crossflow is not used for producers
GEOS_THROW_IF( isProducer() && m_isCrossflowEnabled == 0,
getWrapperDataContext( viewKeyStruct::enableCrossflowString() ) <<
": This option cannot be set to '0' for producers",
InputError );

// 8) Make sure that the initial pressure coefficient is positive
GEOS_THROW_IF( m_initialPressureCoefficient < 0,
getWrapperDataContext( viewKeyStruct::initialPressureCoefficientString() ) <<
Expand Down
108 changes: 54 additions & 54 deletions src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/PerforationFluxKernels.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -103,7 +103,8 @@ class PerforationFluxKernel
CompFlowAccessors const & compFlowAccessors,
MultiFluidAccessors const & multiFluidAccessors,
RelPermAccessors const & relPermAccessors,
bool const disableReservoirToWellFlow ):
bool const isInjector,
bool const isCrossflowEnabled ):
m_resPres( compFlowAccessors.get( fields::flow::pressure {} )),
m_resPhaseVolFrac( compFlowAccessors.get( fields::flow::phaseVolumeFraction {} )),
m_dResPhaseVolFrac( compFlowAccessors.get( fields::flow::dPhaseVolumeFraction {} )),
Expand Down Expand Up @@ -132,7 +133,8 @@ class PerforationFluxKernel
m_compPerfRate( perforationData->getField< fields::well::compPerforationRate >()),
m_dCompPerfRate( perforationData->getField< fields::well::dCompPerforationRate >()),
m_perfStatus( perforationData->getField< fields::perforation::perforationStatus >()),
m_disableReservoirToWellFlow( disableReservoirToWellFlow )
m_isInjector( isInjector ),
m_isCrossflowEnabled( isCrossflowEnabled )
{}

template< typename FUNC = NoOpFunc >
Expand Down Expand Up @@ -254,7 +256,7 @@ class PerforationFluxKernel
// skip the rest of the calculation if the phase is absent
// or if crossflow is disabled for injectors
bool const phaseExists = (m_resPhaseVolFrac[er][esr][ei][ip] > 0);
if( !phaseExists || m_disableReservoirToWellFlow )
if( !phaseExists || (m_isInjector && !m_isCrossflowEnabled) )
{
continue;
}
Expand Down Expand Up @@ -376,8 +378,9 @@ class PerforationFluxKernel
{

// skip the rest of the calculation if the phase is absent
// or if crossflow is disabled for non-injectors aka producers
bool const phaseExists = (m_resPhaseVolFrac[er][esr][ei][ip] > 0);
if( !phaseExists )
if( !phaseExists || (!m_isInjector && !m_isCrossflowEnabled) )
{
continue;
}
Expand Down Expand Up @@ -534,8 +537,9 @@ class PerforationFluxKernel
arrayView3d< real64 > const m_dCompPerfRate_dPres;
arrayView4d< real64 > const m_dCompPerfRate_dComp;
arrayView1d< integer > const m_perfStatus;
bool const m_disableReservoirToWellFlow;

bool const m_isInjector;
bool const m_isCrossflowEnabled;

};

Expand All @@ -550,14 +554,13 @@ class PerforationFluxKernelFactory
* @brief Create a new kernel and launch
* @tparam POLICY the policy used in the RAJA kernel
* @param[in] numComps the number of fluid components
* @param[in] dt time step size
* @param[in] rankOffset the offset of my MPI rank
* @param[in] useTotalMassEquation flag specifying whether to replace one component bal eqn with total mass eqn
* @param[in] dofKey string to get the element degrees of freedom numbers
* @param[in] wellControls object holding well control/constraint information
* @param[in] numPhases the number of phases
* @param[in] flowSolverName the name of the flow solver
* @param[in] perforationData the perforation data
* @param[in] subregion well subregion
* @param[inout] localMatrix the local CRS matrix
* @param[inout] localRhs the local right-hand side vector
* @param[in] elemManager element region manager
* @param[in] isInjector flag indicating whether the well is an injector
* @param[in] isCrossflowEnabled flag indicating whether crossflow is enabled
*/
template< typename POLICY >
static void
Expand All @@ -567,7 +570,8 @@ class PerforationFluxKernelFactory
PerforationData * const perforationData,
ElementSubRegionBase const & subRegion,
ElementRegionManager & elemManager,
integer const disableReservoirToWellFlow )
bool const isInjector,
bool const isCrossflowEnabled )
{
geos::internal::kernelLaunchSelectorCompPhaseSwitch( numComp, numPhases, [&]( auto NC, auto NP )
{
Expand All @@ -580,7 +584,9 @@ class PerforationFluxKernelFactory
typename kernelType::MultiFluidAccessors multiFluidAccessors( elemManager, flowSolverName );
typename kernelType::RelPermAccessors relPermAccessors( elemManager, flowSolverName );

kernelType kernel( perforationData, subRegion, compFlowAccessors, multiFluidAccessors, relPermAccessors, disableReservoirToWellFlow );
kernelType kernel( perforationData, subRegion,
compFlowAccessors, multiFluidAccessors, relPermAccessors,
isInjector, isCrossflowEnabled );
kernelType::template launch< POLICY >( perforationData->size(), kernel );
} );
}
Expand All @@ -601,13 +607,10 @@ class PerforationFluxKernel : public isothermalPerforationFluxKernels::Perforati
public:

using Base = isothermalPerforationFluxKernels::PerforationFluxKernel< NC, NP, IS_THERMAL >;
//using AbstractBase::m_dPhaseVolFrac;
using Base::m_resPhaseCompFrac;
using Base::m_dResCompFrac_dCompDens;
using Base::m_dWellElemCompFrac_dCompDens;

//using AbstractBase::m_dPhaseCompFrac;
//using AbstractBase::m_dCompFrac_dCompDens;
/// Compile time value for the number of components
static constexpr integer numComp = NC;

Expand All @@ -631,10 +634,6 @@ class PerforationFluxKernel : public isothermalPerforationFluxKernels::Perforati
fields::multifluid::phaseEnthalpy,
fields::multifluid::dPhaseEnthalpy >;

//using ThermalConductivityAccessors =
// StencilMaterialAccessors< MultiPhaseThermalConductivityBase,
// fields::thermalconductivity::effectiveConductivity >;

/**
* @brief The type for element-based non-constitutive data parameters.
* Consists entirely of ArrayView's.
Expand All @@ -647,29 +646,30 @@ class PerforationFluxKernel : public isothermalPerforationFluxKernels::Perforati

PerforationFluxKernel ( PerforationData * const perforationData,
ElementSubRegionBase const & subRegion,
MultiFluidBase const & fluid,
MultiFluidBase const & wellFluid,
CompFlowAccessors const & compFlowAccessors,
MultiFluidAccessors const & multiFluidAccessors,
RelPermAccessors const & relPermAccessors,
bool const disableReservoirToWellFlow,
ThermalCompFlowAccessors const & thermalCompFlowAccessors,
ThermalMultiFluidAccessors const & thermalMultiFluidAccessors )
ThermalMultiFluidAccessors const & thermalMultiFluidAccessors,
bool const isInjector,
bool const isCrossflowEnabled )
: Base( perforationData,
subRegion,
compFlowAccessors,
multiFluidAccessors,
relPermAccessors,
disableReservoirToWellFlow ),
m_wellElemPhaseFrac( fluid.phaseFraction() ),
m_dPhaseFrac( fluid.dPhaseFraction() ),
m_wellElemPhaseEnthalpy( fluid.phaseEnthalpy()),
m_dWellElemPhaseEnthalpy( fluid.dPhaseEnthalpy()),
m_energyPerfFlux( perforationData->getField< fields::well::energyPerforationFlux >()),
m_dEnergyPerfFlux( perforationData->getField< fields::well::dEnergyPerforationFlux >()),
isInjector,
isCrossflowEnabled ),
m_wellElemPhaseFrac( wellFluid.phaseFraction() ),
m_dWellElemPhaseFrac( wellFluid.dPhaseFraction() ),
m_wellElemPhaseEnthalpy( wellFluid.phaseEnthalpy()),
m_dWellElemPhaseEnthalpy( wellFluid.dPhaseEnthalpy()),
m_energyPerfFlux( perforationData->getField< fields::well::energyPerforationFlux >() ),
m_dEnergyPerfFlux( perforationData->getField< fields::well::dEnergyPerforationFlux >() ),
m_temp( thermalCompFlowAccessors.get( fields::flow::temperature {} ) ),
m_resPhaseEnthalpy( thermalMultiFluidAccessors.get( fields::multifluid::phaseEnthalpy {} ) ),
m_dResPhaseEnthalpy( thermalMultiFluidAccessors.get( fields::multifluid::dPhaseEnthalpy {} ) )

{}

template< typename FUNC = NoOpFunc >
Expand Down Expand Up @@ -730,32 +730,32 @@ class PerforationFluxKernel : public isothermalPerforationFluxKernels::Perforati
bool const phaseExists = m_wellElemPhaseFrac[iwelem][0][iphase] > 0.0;
if( !phaseExists )
continue;
double pflux = m_wellElemPhaseFrac[iwelem][0][iphase]*flux;
real64 const phaseFlux = m_wellElemPhaseFrac[iwelem][0][iphase] * flux;
real64 const wellelem_enthalpy = m_wellElemPhaseEnthalpy[iwelem][0][iphase];
m_energyPerfFlux[iperf] += pflux * wellelem_enthalpy;
m_energyPerfFlux[iperf] += phaseFlux * wellelem_enthalpy;

// energy equation derivatives WRT res P & T
m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dP] += dFlux[TAG::RES][CP_Deriv::dP] * wellelem_enthalpy;
m_dEnergyPerfFlux[iperf][TAG::RES][CP_Deriv::dT] += dFlux[TAG::RES][CP_Deriv::dT] * wellelem_enthalpy;

m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dP] += dFlux[TAG::WELL][CP_Deriv::dP] * wellelem_enthalpy
+ pflux * m_dWellElemPhaseEnthalpy[iwelem][0][iphase][Deriv::dP]
+ pflux * wellelem_enthalpy * m_dPhaseFrac[iwelem][0][iphase][Deriv::dP];
+ phaseFlux * m_dWellElemPhaseEnthalpy[iwelem][0][iphase][Deriv::dP]
+ phaseFlux * wellelem_enthalpy * m_dWellElemPhaseFrac[iwelem][0][iphase][Deriv::dP];
m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dT] += dFlux[TAG::WELL][CP_Deriv::dT] * wellelem_enthalpy
+ pflux * m_dWellElemPhaseEnthalpy[iwelem][0][iphase][Deriv::dT]
+ pflux * wellelem_enthalpy * m_dPhaseFrac[iwelem][0][iphase][Deriv::dT];
+ phaseFlux * m_dWellElemPhaseEnthalpy[iwelem][0][iphase][Deriv::dT]
+ phaseFlux * wellelem_enthalpy * m_dWellElemPhaseFrac[iwelem][0][iphase][Deriv::dT];

//energy e
real64 dPVF_dC[numComp]{};
applyChainRule( NC,
m_dWellElemCompFrac_dCompDens[iwelem],
m_dPhaseFrac[iwelem][0][iphase],
m_dWellElemPhaseFrac[iwelem][0][iphase],
dPVF_dC,
Deriv::dC );
for( integer ic=0; ic<NC; ic++ )
{
m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dC+ic] += wellelem_enthalpy * dFlux[TAG::WELL][CP_Deriv::dC+ic] * m_wellElemPhaseFrac[iwelem][0][iphase];
m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dC+ic] += wellelem_enthalpy * pflux * dPVF_dC[ ic];
m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dC+ic] += wellelem_enthalpy * dFlux[TAG::WELL][CP_Deriv::dC+ic] * m_wellElemPhaseFrac[iwelem][0][iphase];
m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dC+ic] += wellelem_enthalpy * phaseFlux * dPVF_dC[ic];
}
// energy equation enthalpy derivatives WRT well dens
real64 dProp_dC[numComp]{};
Expand All @@ -767,7 +767,7 @@ class PerforationFluxKernel : public isothermalPerforationFluxKernels::Perforati

for( integer jc = 0; jc < NC; ++jc )
{
m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dC+jc] += pflux * dProp_dC[jc];
m_dEnergyPerfFlux[iperf][TAG::WELL][CP_Deriv::dC+jc] += phaseFlux * dProp_dC[jc];
}
}

Expand Down Expand Up @@ -802,7 +802,7 @@ class PerforationFluxKernel : public isothermalPerforationFluxKernels::Perforati
/// Views on well element properties
/// Element phase fraction
arrayView3d< real64 const, multifluid::USD_PHASE > const m_wellElemPhaseFrac;
arrayView4d< real64 const, multifluid::USD_PHASE_DC > const m_dPhaseFrac;
arrayView4d< real64 const, multifluid::USD_PHASE_DC > const m_dWellElemPhaseFrac;
arrayView3d< real64 const, multifluid::USD_PHASE > const m_wellElemPhaseEnthalpy;
arrayView4d< real64 const, multifluid::USD_PHASE_DC > const m_dWellElemPhaseEnthalpy;

Expand Down Expand Up @@ -831,14 +831,14 @@ class PerforationFluxKernelFactory
* @brief Create a new kernel and launch
* @tparam POLICY the policy used in the RAJA kernel
* @param[in] numComps the number of fluid components
* @param[in] dt time step size
* @param[in] rankOffset the offset of my MPI rank
* @param[in] useTotalMassEquation flag specifying whether to replace one component bal eqn with total mass eqn
* @param[in] dofKey string to get the element degrees of freedom numbers
* @param[in] wellControls object holding well control/constraint information
* @param[in] numPhases the number of phases
* @param[in] flowSolverName the name of the flow solver
* @param[in] perforationData the perforation data
* @param[in] subregion well subregion
* @param[inout] localMatrix the local CRS matrix
* @param[inout] localRhs the local right-hand side vector
* @param[in] fluid the multifluid object
* @param[in] elemManager element region manager
* @param[in] isInjector flag indicating whether the well is an injector
* @param[in] isCrossflowEnabled flag indicating whether crossflow is enabled
*/
template< typename POLICY >
static void
Expand All @@ -849,7 +849,8 @@ class PerforationFluxKernelFactory
ElementSubRegionBase const & subRegion,
MultiFluidBase const & fluid,
ElementRegionManager & elemManager,
integer const disableReservoirToWellFlow )
bool const isInjector,
bool const isCrossflowEnabled )
{
geos::internal::kernelLaunchSelectorCompPhaseSwitch( numComp, numPhases, [&]( auto NC, auto NP )
{
Expand All @@ -865,9 +866,8 @@ class PerforationFluxKernelFactory
typename kernelType::ThermalMultiFluidAccessors thermalMultiFluidAccessors( elemManager, flowSolverName );

kernelType kernel( perforationData, subRegion, fluid, compFlowAccessors, multiFluidAccessors,
relPermAccessors, disableReservoirToWellFlow,
thermalCompFlowAccessors,
thermalMultiFluidAccessors );
relPermAccessors, thermalCompFlowAccessors, thermalMultiFluidAccessors,
isInjector, isCrossflowEnabled );
kernelType::template launch< POLICY >( perforationData->size(), kernel );
} );
}
Expand Down
Loading