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"/> 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 6414693bf22..40eb3528bfe 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,13 @@ 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; + bool usePTDer = false; if( useSurfaceConditions ) { // use surface conditions flashPressure = wellControls.getSurfacePressure(); + flashTemperature = wellControls.getSurfaceTemperature(); } else { @@ -309,10 +313,13 @@ 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]; + usePTDer = true; } } real64 & currentVolRate = @@ -332,12 +339,15 @@ 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, + &usePTDer, &flashPressure, + &flashTemperature, ¤tVolRate, dCurrentVolRate, &iwelemRef, @@ -350,12 +360,14 @@ void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & e if( useSurfaceConditions ) { // we need to compute the surface density - fluidWrapper.update( iwelemRef, 0, flashPressure ); + + fluidWrapper.update( iwelemRef, 0, flashPressure, flashTemperature ); + 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,17 +378,20 @@ void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & e else { real64 const refPres = pres[iwelemRef]; - fluidWrapper.update( iwelemRef, 0, refPres ); + real64 const refTemp = temp[iwelemRef]; + + fluidWrapper.update( iwelemRef, 0, refPres, refTemp ); + } 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 ) { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.cpp index e04bd780afa..4fd4adca812 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()) &&