From 88665f12b19e62a47c1e06c29eacefd12f438c11 Mon Sep 17 00:00:00 2001 From: Tom Byer <149726499+tjb-ltk@users.noreply.github.com> Date: Fri, 20 Feb 2026 09:05:07 -0800 Subject: [PATCH 1/5] call flash with P&T --- .../fluidFlow/wells/SinglePhaseWell.cpp | 32 ++++++++++++++++--- 1 file changed, 27 insertions(+), 5 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp index 6414693bf22..d4999d44d2e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp @@ -266,6 +266,8 @@ void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & e arrayView1d< real64 const > const pres = subRegion.getField< well::pressure >(); + arrayView1d< real64 const > const & temp = subRegion.getField< well::temperature >(); + arrayView1d< real64 const > const & connRate = subRegion.getField< well::connectionRate >(); @@ -282,11 +284,12 @@ void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & e string const wellControlsName = wellControls.getName(); bool const logSurfaceCondition = isLogLevelActive< logInfo::WellControl >( wellControls.getLogLevel()); integer const useSurfaceConditions = wellControls.useSurfaceConditions(); - real64 flashPressure; + real64 flashPressure, flashTemperature; if( useSurfaceConditions ) { // use surface conditions flashPressure = wellControls.getSurfacePressure(); + flashTemperature = wellControls.getSurfaceTemperature(); } else { @@ -309,10 +312,12 @@ void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & e } // use region conditions flashPressure = wellControls.getRegionAveragePressure(); + flashTemperature = wellControls.getRegionAverageTemperature(); if( flashPressure < 0.0 ) { // use segment conditions flashPressure = pres[iwelemRef]; + flashTemperature = temp[iwelemRef]; } } real64 & currentVolRate = @@ -332,12 +337,14 @@ void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & e // bring everything back to host, capture the scalars by reference forAll< serialPolicy >( 1, [fluidWrapper, pres, + temp, connRate, dens, dDens, logSurfaceCondition, &useSurfaceConditions, &flashPressure, + &flashTemperature, ¤tVolRate, dCurrentVolRate, &iwelemRef, @@ -350,12 +357,19 @@ void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & e if( useSurfaceConditions ) { // we need to compute the surface density - fluidWrapper.update( iwelemRef, 0, flashPressure ); + if constexpr ( IS_THERMAL ) + { + fluidWrapper.update( iwelemRef, 0, flashPressure, flashTemperature ); + } + else + { + fluidWrapper.update( iwelemRef, 0, flashPressure ); + } if( logSurfaceCondition ) { - GEOS_LOG_RANK( GEOS_FMT( "{}: surface density computed with P_surface = {} Pa", - wellControlsName, flashPressure ) ); + GEOS_LOG_RANK( GEOS_FMT( "{}: surface density computed with P_surface = {} Pa, T_surface = {} K", + wellControlsName, flashPressure, flashTemperature ) ); } #ifdef GEOS_USE_HIP @@ -366,7 +380,15 @@ void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & e else { real64 const refPres = pres[iwelemRef]; - fluidWrapper.update( iwelemRef, 0, refPres ); + real64 const refTemp = temp[iwelemRef]; + if constexpr ( IS_THERMAL ) + { + fluidWrapper.update( iwelemRef, 0, refPres, refTemp ); + } + else + { + fluidWrapper.update( iwelemRef, 0, refPres, refTemp ); + } } real64 const densInv = 1.0 / dens[iwelemRef][0]; From c051dd7a4744c0f32e809adf175f6c5af46c2325 Mon Sep 17 00:00:00 2001 From: Tom Byer <149726499+tjb-ltk@users.noreply.github.com> Date: Tue, 24 Feb 2026 10:51:38 -0800 Subject: [PATCH 2/5] pr comment requests --- .../fluidFlow/wells/SinglePhaseWell.cpp | 22 +++++-------------- .../fluidFlow/wells/WellControls.cpp | 15 +++++++++++-- 2 files changed, 19 insertions(+), 18 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp index d4999d44d2e..0de3205d8ea 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp @@ -357,14 +357,9 @@ void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & e if( useSurfaceConditions ) { // we need to compute the surface density - if constexpr ( IS_THERMAL ) - { - fluidWrapper.update( iwelemRef, 0, flashPressure, flashTemperature ); - } - else - { - fluidWrapper.update( iwelemRef, 0, flashPressure ); - } + + fluidWrapper.update( iwelemRef, 0, flashPressure, flashTemperature ); + if( logSurfaceCondition ) { @@ -381,14 +376,9 @@ void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & e { real64 const refPres = pres[iwelemRef]; real64 const refTemp = temp[iwelemRef]; - if constexpr ( IS_THERMAL ) - { - fluidWrapper.update( iwelemRef, 0, refPres, refTemp ); - } - else - { - fluidWrapper.update( iwelemRef, 0, refPres, refTemp ); - } + + fluidWrapper.update( iwelemRef, 0, refPres, refTemp ); + } real64 const densInv = 1.0 / dens[iwelemRef][0]; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp index 684306cce46..0a5535a46b5 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp @@ -41,8 +41,8 @@ WellControls::WellControls( string const & name, Group * const parent ) m_targetPhaseRate( 0.0 ), m_targetMassRate( 0.0 ), m_useSurfaceConditions( 0 ), - m_surfacePres( 0.0 ), - m_surfaceTemp( 0.0 ), + m_surfacePres( -1.0 ), + m_surfaceTemp( -1.0 ), m_isCrossflowEnabled( 1 ), m_initialPressureCoefficient( 0.1 ), m_rateSign( -1.0 ), @@ -305,6 +305,17 @@ void WellControls::postInputInitialization() "The flag to select surface/reservoir conditions must be equal to 0 or 1", InputError, getWrapperDataContext( viewKeyStruct::useSurfaceConditionsString() ) ); + if( m_useSurfaceConditions ) + { + // 3a) check the flag for surface condition P&T + GEOS_THROW_IF( m_surfacePres <= 0.0 , + "Surface conditions set to 1 but surface pressure is not set", + InputError, getWrapperDataContext( viewKeyStruct::useSurfaceConditionsString() ) ); + GEOS_THROW_IF( m_surfaceTemp <= 0.0 , + "Surface conditions set to 1 but surface temperature is not set", + InputError, getWrapperDataContext( viewKeyStruct::useSurfaceConditionsString() ) ); + + } // 4) check that at least one rate constraint has been defined GEOS_THROW_IF( ((m_targetPhaseRate <= 0.0 && m_targetPhaseRateTableName.empty()) && (m_targetMassRate <= 0.0 && m_targetMassRateTableName.empty()) && From 7e469c75f6e6dff7c20fbf91da09469b22b6ebd8 Mon Sep 17 00:00:00 2001 From: Tom Byer <149726499+tjb-ltk@users.noreply.github.com> Date: Tue, 24 Feb 2026 12:04:31 -0800 Subject: [PATCH 3/5] uncrust --- .../physicsSolvers/fluidFlow/wells/WellControls.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp index 0a5535a46b5..120a28352bc 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp @@ -308,10 +308,10 @@ void WellControls::postInputInitialization() if( m_useSurfaceConditions ) { // 3a) check the flag for surface condition P&T - GEOS_THROW_IF( m_surfacePres <= 0.0 , + GEOS_THROW_IF( m_surfacePres <= 0.0, "Surface conditions set to 1 but surface pressure is not set", InputError, getWrapperDataContext( viewKeyStruct::useSurfaceConditionsString() ) ); - GEOS_THROW_IF( m_surfaceTemp <= 0.0 , + GEOS_THROW_IF( m_surfaceTemp <= 0.0, "Surface conditions set to 1 but surface temperature is not set", InputError, getWrapperDataContext( viewKeyStruct::useSurfaceConditionsString() ) ); From f2dd4c070f4c7297e309445b31bd495a168ae00d Mon Sep 17 00:00:00 2001 From: Tom Byer <149726499+tjb-ltk@users.noreply.github.com> Date: Wed, 25 Feb 2026 04:42:30 -0800 Subject: [PATCH 4/5] add temp to sep def for wells --- .../poromechanics/PoroElastic_staircase_singlephase_3d_fim.xml | 1 + .../PoroElastic_staircase_singlephase_3d_sequential.xml | 1 + .../singlePhaseWell/staircase_single_phase_wells_hybrid_3d.xml | 1 + 3 files changed, 3 insertions(+) diff --git a/inputFiles/poromechanics/PoroElastic_staircase_singlephase_3d_fim.xml b/inputFiles/poromechanics/PoroElastic_staircase_singlephase_3d_fim.xml index 95f2090ada1..395292ba46d 100644 --- a/inputFiles/poromechanics/PoroElastic_staircase_singlephase_3d_fim.xml +++ b/inputFiles/poromechanics/PoroElastic_staircase_singlephase_3d_fim.xml @@ -63,6 +63,7 @@ targetBHP="5e7" targetTotalRateTableName="injectorTotalRateTable" useSurfaceConditions="1" + surfaceTemperature="300" surfacePressure="101325"/> diff --git a/inputFiles/poromechanics/PoroElastic_staircase_singlephase_3d_sequential.xml b/inputFiles/poromechanics/PoroElastic_staircase_singlephase_3d_sequential.xml index 0c952c6934f..af1007c6c6a 100755 --- a/inputFiles/poromechanics/PoroElastic_staircase_singlephase_3d_sequential.xml +++ b/inputFiles/poromechanics/PoroElastic_staircase_singlephase_3d_sequential.xml @@ -72,6 +72,7 @@ targetBHP="5e7" targetTotalRateTableName="injectorTotalRateTable" useSurfaceConditions="1" + surfaceTemperature="300" surfacePressure="101325"/> diff --git a/inputFiles/singlePhaseWell/staircase_single_phase_wells_hybrid_3d.xml b/inputFiles/singlePhaseWell/staircase_single_phase_wells_hybrid_3d.xml index 2298a32e483..52f243e506d 100644 --- a/inputFiles/singlePhaseWell/staircase_single_phase_wells_hybrid_3d.xml +++ b/inputFiles/singlePhaseWell/staircase_single_phase_wells_hybrid_3d.xml @@ -38,6 +38,7 @@ control="totalVolRate" useSurfaceConditions="1" surfacePressure="101325" + surfaceTemperature="300" referenceElevation="12" targetBHP="1e7" targetTotalRate="1e-6"/> From d31f70c0900ac20811e216e4f89ebab5a20ecab8 Mon Sep 17 00:00:00 2001 From: Tom Byer <149726499+tjb-ltk@users.noreply.github.com> Date: Tue, 10 Mar 2026 15:04:15 -0700 Subject: [PATCH 5/5] resvol constraint der fix --- .../fluidFlow/wells/CompositionalMultiphaseWell.cpp | 11 +++++++---- .../fluidFlow/wells/SinglePhaseWell.cpp | 7 +++++-- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index 275b0451d82..fb88119b589 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -709,6 +709,7 @@ void CompositionalMultiphaseWell::updateVolRatesForConstraint( ElementRegionMana integer const useSurfaceConditions = wellControls.useSurfaceConditions(); real64 flashPressure; real64 flashTemperature; + bool usePTDer = false; if( useSurfaceConditions ) { // use surface conditions @@ -741,6 +742,7 @@ void CompositionalMultiphaseWell::updateVolRatesForConstraint( ElementRegionMana // region name not set, use segment conditions flashPressure = pres[iwelemRef]; flashTemperature = temp[iwelemRef]; + usePTDer = true; } else { @@ -790,6 +792,7 @@ void CompositionalMultiphaseWell::updateVolRatesForConstraint( ElementRegionMana dPhaseFrac, logSurfaceCondition, &useSurfaceConditions, + &usePTDer, &flashPressure, &flashTemperature, ¤tTotalVolRate, @@ -849,10 +852,10 @@ void CompositionalMultiphaseWell::updateVolRatesForConstraint( ElementRegionMana currentTotalVolRate = currentTotalRate * totalDensInv; // Compute derivatives dP dT real64 const dTotalDensInv_dPres = -dTotalDens[iwelemRef][0][Deriv::dP] * totalDensInv * totalDensInv; - dCurrentTotalVolRate[COFFSET_WJ::dP] = ( useSurfaceConditions == 0 ) * currentTotalRate * dTotalDensInv_dPres; + dCurrentTotalVolRate[COFFSET_WJ::dP] = ( usePTDer == 1 ) * currentTotalRate * dTotalDensInv_dPres; if constexpr ( IS_THERMAL ) { - dCurrentTotalVolRate[COFFSET_WJ::dT] = ( useSurfaceConditions == 0 ) * currentTotalRate * -dTotalDens[iwelemRef][0][Deriv::dT] * totalDensInv * totalDensInv; + dCurrentTotalVolRate[COFFSET_WJ::dT] = ( usePTDer == 1 ) * currentTotalRate * -dTotalDens[iwelemRef][0][Deriv::dT] * totalDensInv * totalDensInv; } if( logSurfaceCondition && useSurfaceConditions ) @@ -888,13 +891,13 @@ void CompositionalMultiphaseWell::updateVolRatesForConstraint( ElementRegionMana // Step 3.2: divide the total mass/molar rate by the (phase density * phase fraction) to get the phase volumetric rate currentPhaseVolRate[ip] = currentTotalRate * phaseFracTimesPhaseDensInv; - dCurrentPhaseVolRate[ip][COFFSET_WJ::dP] = ( useSurfaceConditions == 0 ) * currentTotalRate * dPhaseFracTimesPhaseDensInv_dPres; + dCurrentPhaseVolRate[ip][COFFSET_WJ::dP] = ( usePTDer == 1 ) * currentTotalRate * dPhaseFracTimesPhaseDensInv_dPres; dCurrentPhaseVolRate[ip][COFFSET_WJ::dQ] = phaseFracTimesPhaseDensInv; if constexpr (IS_THERMAL ) { real64 const dPhaseFracTimesPhaseDensInv_dTemp = dPhaseFrac[iwelemRef][0][ip][Deriv::dT] * phaseDensInv - dPhaseDens[iwelemRef][0][ip][Deriv::dT] * phaseFracTimesPhaseDensInv * phaseDensInv; - dCurrentPhaseVolRate[ip][COFFSET_WJ::dT] = ( useSurfaceConditions == 0 ) * currentTotalRate * dPhaseFracTimesPhaseDensInv_dTemp; + dCurrentPhaseVolRate[ip][COFFSET_WJ::dT] = ( usePTDer == 1 ) * currentTotalRate * dPhaseFracTimesPhaseDensInv_dTemp; } for( integer ic = 0; ic < numComp; ++ic ) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp index 0de3205d8ea..40eb3528bfe 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp @@ -285,6 +285,7 @@ void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & e bool const logSurfaceCondition = isLogLevelActive< logInfo::WellControl >( wellControls.getLogLevel()); integer const useSurfaceConditions = wellControls.useSurfaceConditions(); real64 flashPressure, flashTemperature; + bool usePTDer = false; if( useSurfaceConditions ) { // use surface conditions @@ -318,6 +319,7 @@ void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & e // use segment conditions flashPressure = pres[iwelemRef]; flashTemperature = temp[iwelemRef]; + usePTDer = true; } } real64 & currentVolRate = @@ -343,6 +345,7 @@ void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & e dDens, logSurfaceCondition, &useSurfaceConditions, + &usePTDer, &flashPressure, &flashTemperature, ¤tVolRate, @@ -384,11 +387,11 @@ void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & e real64 const densInv = 1.0 / dens[iwelemRef][0]; currentVolRate = connRate[iwelemRef] * densInv; - dCurrentVolRate[COFFSET_WJ::dP] = -( useSurfaceConditions == 0 ) * dDens[iwelemRef][0][DerivOffset::dP] * currentVolRate * densInv; + dCurrentVolRate[COFFSET_WJ::dP] = -( usePTDer == 1 ) * dDens[iwelemRef][0][DerivOffset::dP] * currentVolRate * densInv; dCurrentVolRate[COFFSET_WJ::dQ] = densInv; if constexpr ( IS_THERMAL ) { - dCurrentVolRate[COFFSET_WJ::dT] = -( useSurfaceConditions == 0 ) * dDens[iwelemRef][0][DerivOffset::dT] * currentVolRate * densInv; + dCurrentVolRate[COFFSET_WJ::dT] = -( usePTDer == 1 ) * dDens[iwelemRef][0][DerivOffset::dT] * currentVolRate * densInv; } if( logSurfaceCondition && useSurfaceConditions ) {