diff --git a/inputFiles/singlePhaseWell/thermal_singlePhase_well.xml b/inputFiles/singlePhaseWell/thermal_singlePhase_well.xml new file mode 100755 index 00000000000..78d27f5283c --- /dev/null +++ b/inputFiles/singlePhaseWell/thermal_singlePhase_well.xml @@ -0,0 +1,266 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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() diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndSinglePhaseWellKernels.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndSinglePhaseWellKernels.hpp index b16360f690a..37f006aff38 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 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 ); @@ -381,22 +380,22 @@ 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; // 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..fd260fd63e8 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 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]; }