From f5ac0a64814686af239230fe40d883a714935159 Mon Sep 17 00:00:00 2001 From: npillardou Date: Fri, 13 Feb 2026 23:15:02 +0100 Subject: [PATCH 1/6] Add fix for top element well --- ...pledReservoirAndSinglePhaseWellKernels.hpp | 21 ++++++++--------- .../CoupledReservoirAndWellKernels.hpp | 23 +++++++++---------- 2 files changed, 21 insertions(+), 23 deletions(-) diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndSinglePhaseWellKernels.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndSinglePhaseWellKernels.hpp index b16360f690a..a573a62266b 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndSinglePhaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndSinglePhaseWellKernels.hpp @@ -366,14 +366,13 @@ class ThermalSinglePhaseFluxKernel : public IsothermalSinglePhaseFluxKernel< IS_ stackArray1d< globalIndex, 2*resNumDOF > & dofColIndices, localIndex const iwelem ) { - // No energy equation if top element and Injector - // Top element defined by global index == 0 - // Assumption is global index == 0 is top segment with fixed temp BC - if( !m_isProducer ) - { - if( m_globalWellElementIndex[iwelem] == 0 ) - return; - } + // For injector top element (global index == 0), the well energy equation + // is replaced by a Dirichlet BC (T = T_inj), so we must NOT assemble + // the well-side energy flux. However, we MUST still assemble the + // reservoir-side energy flux so the reservoir cell receives the correct + // enthalpy from the injected mass. + bool const isTopInjectorElement = !m_isProducer && m_globalWellElementIndex[iwelem] == 0; + // local working variables and arrays stackArray1d< localIndex, 2 > eqnRowIndices( 2 ); @@ -387,16 +386,16 @@ class ThermalSinglePhaseFluxKernel : public IsothermalSinglePhaseFluxKernel< IS_ // populate local flux vector and derivatives localPerf[TAG::RES ] = m_dt * m_energyPerfFlux[iperf]; - localPerf[TAG::WELL ] = -m_dt * m_energyPerfFlux[iperf]; + localPerf[TAG::WELL ] = isTopInjectorElement ? 0.0 : -m_dt * m_energyPerfFlux[iperf]; for( integer ke = 0; ke < 2; ++ke ) { localIndex localDofIndexPres = ke * resNumDOF; localPerfJacobian[TAG::RES ][localDofIndexPres] = m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dP]; - localPerfJacobian[TAG::WELL ][localDofIndexPres] = -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dP]; + localPerfJacobian[TAG::WELL ][localDofIndexPres] = isTopInjectorElement ? 0.0 : -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dP]; localPerfJacobian[TAG::RES ][localDofIndexPres+1] = m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dT]; - localPerfJacobian[TAG::WELL][localDofIndexPres+1] = -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dT]; + localPerfJacobian[TAG::WELL][localDofIndexPres+1] = isTopInjectorElement ? 0.0 : -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dT]; } diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellKernels.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellKernels.hpp index 5d106df3156..16307ca4330 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellKernels.hpp @@ -432,14 +432,13 @@ class ThermalCompositionalMultiPhaseFluxKernel : public IsothermalCompositionalM stackArray1d< globalIndex, 2*resNumDOF > & dofColIndices, localIndex const iwelem ) { - // No energy equation if top element and Injector - // Top element defined by global index == 0 - // Assumption is global index == 0 is top segment with fixed temp BC - if( !m_isProducer ) - { - if( m_globalWellElementIndex[iwelem] == 0 ) - return; - } + // For injector top element (global index == 0), the well energy equation + // is replaced by a Dirichlet BC (T = T_inj), so we must NOT assemble + // the well-side energy flux. However, we MUST still assemble the + // reservoir-side energy flux so the reservoir cell receives the correct + // enthalpy from the injected mass. + bool const isTopInjectorElement = !m_isProducer && m_globalWellElementIndex[iwelem] == 0; + // local working variables and arrays stackArray1d< localIndex, 2* numComp > eqnRowIndices( 2 ); @@ -453,23 +452,23 @@ class ThermalCompositionalMultiPhaseFluxKernel : public IsothermalCompositionalM // populate local flux vector and derivatives localPerf[TAG::RES ] = m_dt * m_energyPerfFlux[iperf]; - localPerf[TAG::WELL ] = -m_dt * m_energyPerfFlux[iperf]; + localPerf[TAG::WELL ] = isTopInjectorElement ? 0.0 : -m_dt * m_energyPerfFlux[iperf]; for( integer ke = 0; ke < 2; ++ke ) { localIndex localDofIndexPres = ke * resNumDOF; localPerfJacobian[TAG::RES ][localDofIndexPres] = m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dP]; - localPerfJacobian[TAG::WELL ][localDofIndexPres] = -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dP]; + localPerfJacobian[TAG::WELL ][localDofIndexPres] = isTopInjectorElement ? 0.0 : -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dP]; // populate local flux vector and derivatives for( integer ic = 0; ic < numComp; ++ic ) { localIndex const localDofIndexComp = localDofIndexPres + ic + 1; localPerfJacobian[TAG::RES ][localDofIndexComp] = m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dC+ic]; - localPerfJacobian[TAG::WELL][localDofIndexComp] = -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dC+ic]; + localPerfJacobian[TAG::WELL][localDofIndexComp] = isTopInjectorElement ? 0.0 : -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dC+ic]; } localPerfJacobian[TAG::RES ][localDofIndexPres+NC+1] = m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dT]; - localPerfJacobian[TAG::WELL][localDofIndexPres+NC+1] = -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dT]; + localPerfJacobian[TAG::WELL][localDofIndexPres+NC+1] = isTopInjectorElement ? 0.0 : -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dT]; } From 9d60588148eead7bb68d6fa89441817b2bb3b703 Mon Sep 17 00:00:00 2001 From: npillardou Date: Mon, 16 Feb 2026 09:30:37 +0100 Subject: [PATCH 2/6] add comments --- .../CoupledReservoirAndSinglePhaseWellKernels.hpp | 6 +++--- .../multiphysics/CoupledReservoirAndWellKernels.hpp | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndSinglePhaseWellKernels.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndSinglePhaseWellKernels.hpp index a573a62266b..37f006aff38 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndSinglePhaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndSinglePhaseWellKernels.hpp @@ -367,8 +367,8 @@ class ThermalSinglePhaseFluxKernel : public IsothermalSinglePhaseFluxKernel< IS_ localIndex const iwelem ) { // For injector top element (global index == 0), the well energy equation - // is replaced by a Dirichlet BC (T = T_inj), so we must NOT assemble - // the well-side energy flux. However, we MUST still assemble the + // is replaced by a Dirichlet BC (T = T_inj), so we must not assemble + // the well-side energy flux. However, we still assemble the // reservoir-side energy flux so the reservoir cell receives the correct // enthalpy from the injected mass. bool const isTopInjectorElement = !m_isProducer && m_globalWellElementIndex[iwelem] == 0; @@ -380,7 +380,7 @@ class ThermalSinglePhaseFluxKernel : public IsothermalSinglePhaseFluxKernel< IS_ stackArray2d< real64, 2*2 * resNumDOF > localPerfJacobian( 2, 2 * resNumDOF ); - // equantion offsets - note res and well have different equation lineups + // equation offsets - note res and well have different equation lineups eqnRowIndices[TAG::RES ] = LvArray::integerConversion< localIndex >( resOffset - m_rankOffset ) + 1; eqnRowIndices[TAG::WELL ] = LvArray::integerConversion< localIndex >( wellElemOffset - m_rankOffset ) + WJ_ROFFSET::ENERGYBAL; diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellKernels.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellKernels.hpp index 16307ca4330..fd260fd63e8 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellKernels.hpp @@ -433,8 +433,8 @@ class ThermalCompositionalMultiPhaseFluxKernel : public IsothermalCompositionalM localIndex const iwelem ) { // For injector top element (global index == 0), the well energy equation - // is replaced by a Dirichlet BC (T = T_inj), so we must NOT assemble - // the well-side energy flux. However, we MUST still assemble the + // is replaced by a Dirichlet BC (T = T_inj), so we must not assemble + // the well-side energy flux. However, we still assemble the // reservoir-side energy flux so the reservoir cell receives the correct // enthalpy from the injected mass. bool const isTopInjectorElement = !m_isProducer && m_globalWellElementIndex[iwelem] == 0; From 01d3764db1e99959d8f9aa4917439e54ce707959 Mon Sep 17 00:00:00 2001 From: npillardou Date: Wed, 18 Feb 2026 14:44:38 +0100 Subject: [PATCH 3/6] Add test case for thermal single phase well --- .../thermal_singlePhase_well.xml | 262 ++++++++++++++++++ 1 file changed, 262 insertions(+) create mode 100755 inputFiles/singlePhaseWell/thermal_singlePhase_well.xml diff --git a/inputFiles/singlePhaseWell/thermal_singlePhase_well.xml b/inputFiles/singlePhaseWell/thermal_singlePhase_well.xml new file mode 100755 index 00000000000..e47965af82b --- /dev/null +++ b/inputFiles/singlePhaseWell/thermal_singlePhase_well.xml @@ -0,0 +1,262 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From 3d8b2600b7a1b284a6a57b4e6244b8dbecf62b96 Mon Sep 17 00:00:00 2001 From: npillardou Date: Wed, 18 Feb 2026 17:24:40 +0100 Subject: [PATCH 4/6] Update xml file --- .../thermal_singlePhase_well.xml | 32 ++++++++++--------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/inputFiles/singlePhaseWell/thermal_singlePhase_well.xml b/inputFiles/singlePhaseWell/thermal_singlePhase_well.xml index e47965af82b..7e9b9b09864 100755 --- a/inputFiles/singlePhaseWell/thermal_singlePhase_well.xml +++ b/inputFiles/singlePhaseWell/thermal_singlePhase_well.xml @@ -10,6 +10,7 @@ flowSolverName="SinglePhaseFlow" wellSolverName="SinglePhaseWell" logLevel="1" + initialDt="1" targetRegions="{ Reservoir, well }"> + injectionStream="{1.0}"/> @@ -104,23 +105,21 @@ polylineSegmentConn="{ {0,1}, {1,2}, {2,3}, {3,4}, {4,5} }"> + distanceFromHead="5" /> + distanceFromHead="15" /> + distanceFromHead="25" /> + distanceFromHead="35" /> - + distanceFromHead="45" /> + @@ -231,7 +230,7 @@ @@ -239,7 +238,7 @@ - + + + - + From 519c304caca2957155d313650bdf5b397d868a11 Mon Sep 17 00:00:00 2001 From: npillardou Date: Thu, 19 Feb 2026 15:04:17 +0100 Subject: [PATCH 5/6] Add assertion when no segment above top perfs --- .../mesh/generators/WellGeneratorBase.cpp | 20 ++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/src/coreComponents/mesh/generators/WellGeneratorBase.cpp b/src/coreComponents/mesh/generators/WellGeneratorBase.cpp index f55ba78049d..6e7fbe1ca01 100644 --- a/src/coreComponents/mesh/generators/WellGeneratorBase.cpp +++ b/src/coreComponents/mesh/generators/WellGeneratorBase.cpp @@ -458,6 +458,24 @@ void WellGeneratorBase::checkPerforationLocationsValidity() mergePerforations( elemToPerfMap ); } + // check that the top well element (wellhead) does not have a perforation + // The top element carries the well control equation (BHP or rate constraint). + // Having a perforation on the same element creates a strong coupling between the control constraint + // and the perforation flux, which degrades Newton convergence (both isothermal and thermal). + // The top segment should be reserved for the well boundary condition constraint only. + for( globalIndex iwelem = 0; iwelem < m_numElems; ++iwelem ) + { + GEOS_THROW_IF( m_nextElemId[iwelem] < 0 && elemToPerfMap[iwelem].size() > 0, + "Well '" << getName() << "': a perforation is placed on the top well element (wellhead). " + << "This is not allowed because the top element is reserved for the well control constraint " + << "(BHP or rate) and must not have any perforation. \n\n" + << "To fix this, extend the well polyline upward by adding a segment above the first perforation. " + << "For example, add a node above the current well head in \"polylineNodeCoords\" and a corresponding " + << "entry in \"polylineSegmentConn\". " + << "This top segment will act as a constraint-only segment with no perforation.", + InputError ); + } + for( globalIndex iwelem = 0; iwelem < m_numElems; ++iwelem ) { // check that there is always a perforation in the last well element (otherwise, the problem is not well posed) @@ -504,7 +522,7 @@ void WellGeneratorBase::mergePerforations( array1d< array1d< localIndex > > cons continue; } - GEOS_LOG_RANK_0( "\n \nThe GEOSX wells currently have the following limitation in parallel: \n" + GEOS_LOG_RANK_0( "\n \nThe GEOS wells currently have the following limitation in parallel: \n" << "We cannot allow an element of the well mesh to have two or more perforations associated with it. \n" << "So, in the present simulation, perforation #" << elemToPerfMap[iwelem][ip] << " of well " << getName() From eedc3192593c9b2849a89d2024e973974be600a9 Mon Sep 17 00:00:00 2001 From: npillardou Date: Thu, 19 Feb 2026 16:59:22 +0100 Subject: [PATCH 6/6] Update test case with segments above perf --- .../singlePhaseWell/thermal_singlePhase_well.xml | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/inputFiles/singlePhaseWell/thermal_singlePhase_well.xml b/inputFiles/singlePhaseWell/thermal_singlePhase_well.xml index 7e9b9b09864..78d27f5283c 100755 --- a/inputFiles/singlePhaseWell/thermal_singlePhase_well.xml +++ b/inputFiles/singlePhaseWell/thermal_singlePhase_well.xml @@ -96,28 +96,29 @@ minElementLength="0.0001" minSegmentLength="0.001" logLevel="2" - polylineNodeCoords="{ {495,245,-1600}, + polylineNodeCoords="{ {495,245,-1500}, + {495,245,-1600}, {495,245,-1610}, {495,245,-1620}, {495,245,-1630}, {495,245,-1640}, {495,245,-1650}}" - polylineSegmentConn="{ {0,1}, {1,2}, {2,3}, {3,4}, {4,5} }"> + polylineSegmentConn="{ {0,1}, {1,2}, {2,3}, {3,4}, {4,5}, {5,6} }"> + distanceFromHead="105" /> + distanceFromHead="115" /> + distanceFromHead="125" /> + distanceFromHead="135" /> + distanceFromHead="145" /> @@ -262,3 +263,4 @@ +