Fix thermal well perforation energy flux

inputFiles/singlePhaseWell/thermal_singlePhase_well.xml

@@ -0,0 +1,266 @@
+<?xml version="1.0" ?>
+
+<Problem>
+
+<!-- SPHINX_FIELD_CASE_SOLVER -->
+  <Solvers>
+
+    <SinglePhaseReservoir
+      name="coupledFlowAndWells"
+      flowSolverName="SinglePhaseFlow"
+      wellSolverName="SinglePhaseWell"
+      logLevel="1"
+      initialDt="1"
+      targetRegions="{ Reservoir, well }">
+      <NonlinearSolverParameters
+          newtonTol="1.0e-3"
+          newtonMaxIter="10"
+          maxTimeStepCuts="10"
+          maxSubSteps="100"
+          lineSearchAction="None" />
+      <LinearSolverParameters
+          solverType="fgmres"
+          preconditionerType="mgr"
+          krylovTol="1.0e-4"
+          krylovMaxIter="250"
+          krylovWeakestTol="0.01"
+          logLevel="1" />
+    </SinglePhaseReservoir>
+
+    <SinglePhaseFVM
+      name="SinglePhaseFlow"
+      discretization="TPFA"
+      targetRegions="{ Reservoir }"
+      isThermal="1"
+      temperature="355"
+      logLevel="1">
+      <NonlinearSolverParameters
+        newtonTol="1.0e-6"
+        newtonMaxIter="8"
+        lineSearchAction="None"
+        logLevel="1"/>
+      <LinearSolverParameters
+        solverType="fgmres"
+        preconditionerType="amg"
+        krylovTol="1.0e-10"
+        krylovMaxIter="250"
+        logLevel="1"/>
+    </SinglePhaseFVM>
+
+    <SinglePhaseWell
+      name="SinglePhaseWell"
+      targetRegions="{ well }"
+      isThermal="1"
+      initialDt="86400"
+      writeCSV="1"
+      logLevel="1">
+      <WellControls
+        name="wellControls1"
+        logLevel="1"      
+        type="injector"
+        control="totalVolRate"
+        enableCrossflow="0"
+        referenceElevation="-1650"
+        targetBHP="5e7"
+        useSurfaceConditions="1"
+        surfacePressure="101325"
+        surfaceTemperature="285.15"
+        targetTotalRateTableName="injectorTotalRateTable"     
+        injectionTemperature="293.15"
+        injectionStream="{1.0}"/>
+
+    </SinglePhaseWell>
+
+  </Solvers>
+<!-- SPHINX_FIELD_CASE_SOLVER_END -->
+
+<!-- SPHINX_FIELD_CASE_MESH -->
+<Mesh>
+  <InternalMesh
+    name="mesh1"
+    elementTypes="{ C3D8 }"
+    xCoords="{ 0, 1000 }"
+    yCoords="{ 0, 500 }"
+    zCoords="{ -2000,-1700,-1500, -1000 }"
+    nx="{ 100 }"
+    ny="{ 50 }"
+    nz="{ 10,20,10 }"
+    cellBlockNames="{ cellBlock }">
+
+    <InternalWell
+      name="well"
+      wellRegionName="well"
+      wellControlsName="wellControls1"
+      radius="0.1"
+      numElementsPerSegment="1"
+      minElementLength="0.0001"
+      minSegmentLength="0.001"
+      logLevel="2"
+      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}, {5,6} }">
+      <Perforation
+          name="Perf_1"
+          distanceFromHead="105" />
+      <Perforation
+          name="Perf_2"
+          distanceFromHead="115" />
+      <Perforation
+          name="Perf_3"
+          distanceFromHead="125" />
+      <Perforation
+          name="Perf_4"
+          distanceFromHead="135" />
+      <Perforation
+          name="Perf_5"
+          distanceFromHead="145" />
+    </InternalWell>
+
+  </InternalMesh>
+</Mesh>
+<!-- SPHINX_FIELD_CASE_MESH_END -->
+
+<!-- SPHINX_FIELD_CASE_REGION -->
+  <ElementRegions>
+    <CellElementRegion
+      name="Reservoir"
+      cellBlocks="{ * }"
+      materialList="{ water, rock, thermalCond }"/>
+
+    <WellElementRegion
+        meshBody="mesh1"
+        name="well"
+        materialList="{ water }" />
+
+  </ElementRegions>
+<!-- SPHINX_FIELD_CASE_REGION_END -->
+
+<!-- SPHINX_FIELD_CASE_CONSTITUTIVE -->
+  <Constitutive>
+
+    <ThermalCompressibleSinglePhaseFluid
+        name="water"
+        defaultDensity="1000"
+        defaultViscosity="1.0e-3"
+        referencePressure="1e5"
+        referenceDensity="1000"
+        referenceViscosity="1.0E-3"
+        referenceTemperature="277.15"
+        compressibility="4.5e-10"
+        viscosibility="0.0"
+        thermalExpansionCoeff="1e-5"
+        specificHeatCapacity="4186"
+        referenceInternalEnergy="1"/>
+
+    <SinglePhaseThermalConductivity
+      name="thermalCond"
+      defaultThermalConductivityComponents="{ 2, 2, 2 }"/>
+
+    <CompressibleSolidConstantPermeability
+      name="rock"
+      solidModelName="nullSolid"
+      porosityModelName="rockPorosity"
+      permeabilityModelName="rockPerm"
+      solidInternalEnergyModelName="rockEnergy"/>
+
+    <NullModel
+      name="nullSolid"/>
+
+    <PressurePorosity
+      name="rockPorosity"
+      defaultReferencePorosity="0.05"
+      referencePressure="10e7"
+      compressibility="4.5e-10"/>
+
+    <ConstantPermeability
+      name="rockPerm"
+      permeabilityComponents="{ 1.0e-13, 1.0e-13, 1.0e-15 }"/>
+
+  <SolidInternalEnergy
+      name="rockEnergy"
+      referenceVolumetricHeatCapacity="2.2e6"
+      referenceTemperature="0.0"
+      referenceInternalEnergy="0.0" />
+
+  </Constitutive>
+<!-- SPHINX_FIELD_CASE_CONSTITUTIVE_END -->
+
+<!-- SPHINX_FIELD_CASE_NUMERICAL -->
+  <NumericalMethods>
+    <FiniteVolume>
+      <TwoPointFluxApproximation
+        name="TPFA"/>
+    </FiniteVolume>
+  </NumericalMethods>
+<!-- SPHINX_FIELD_CASE_NUMERICAL_END -->
+
+<!-- SPHINX_FIELD_CASE_FIELD -->
+  <FieldSpecifications>
+    <HydrostaticEquilibrium
+        name="INITIALISATION"
+        datumElevation="-1650"
+        datumPressure="1.65e7"
+        objectPath="ElementRegions"
+        temperatureVsElevationTableName="initTempTable"/>
+
+  </FieldSpecifications>
+<!-- SPHINX_FIELD_CASE_FIELD_END -->
+
+<!-- SPHINX_FIELD_CASE_OUTPUT -->
+  <Outputs>
+    <VTK
+      name="vtkOutput_well1P"
+      fieldNames="{initialPressure, initialTemperature}" />
+  </Outputs>
+<!-- SPHINX_FIELD_CASE_OUTPUT_END -->
+
+<!-- SPHINX_FIELD_CASE_TFUNC -->
+  <Functions>
+    <!-- Geothermal gradient -->
+    <TableFunction
+      name="initTempTable"
+      coordinates="{ -2000, 0.0 }"
+      values="{ 337.15, 277.15 }"
+      interpolation="linear"/>
+
+    <TableFunction
+      name="injectorTotalRateTable"
+      inputVarNames="{ time }"
+      coordinates="{  0, 10e6, 100.0e6 }"
+      values="{ 0, 1e-2, 1e-2 }" 
+      interpolation="lower"/>
+
+  </Functions>
+<!-- SPHINX_FIELD_CASE_TFUNC_END -->
+
+<!-- SPHINX_FIELD_CASE_EVENTS -->
+<Events minTime="0" maxTime="101.0e6">
+
+  <SoloEvent
+    name="Initial_output"
+    targetTime="0"
+    targetExactTimestep="1"
+    target="/Outputs/vtkOutput_well1P"/>
+
+  <PeriodicEvent name="outputs"
+                 beginTime="10.0e6"
+                 timeFrequency="10.0e6"
+                 targetExactTimestep="1"
+                 target="/Outputs/vtkOutput_well1P" />
+
+  <PeriodicEvent name="solverApplications"
+                 maxEventDt="10.0e6"
+                 target="/Solvers/coupledFlowAndWells" />
+
+
+
+</Events>
+<!-- SPHINX_FIELD_CASE_EVENTS_END -->
+
+</Problem>
+

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()

src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndSinglePhaseWellKernels.hpp

@@ -366,37 +366,36 @@ 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 );
 
       stackArray1d< real64, 2 > localPerf( 2 );
       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];
       }
 
 

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];
       }