feat: well_manager + well estimator by tjb-ltk · Pull Request #3972 · GEOS-DEV/GEOS

tjb-ltk · 2026-02-13T15:09:34Z

Replaces PR #3735 and extends PR #3971

Overview

This PR introduces

An improved well schema layout and code refactor better suited for well modeling. The previous implementation primarily targeted Jacobian generation requirements for the coupled reservoir and well system.
The well estimator is used to select the active well constraint by solving the well system assuming fixed reservoir conditions and selecting the constraint with the highest or lowest well flowing pressure.
- the estimator is applied at beginning of Jacobian assembly of the coupled system
- frequency determined via input setting
- estimator is a data member of WellControl and a tailed version of PhysicsSolver for wells
- logic underpins future capability to replace the segmented well model with other formulations

This is a breaking change, detailed migration instructions and migration script are posted at #4081

…broken

…aint class hierachy , 3) Remove all old constraint handling methods/variables, 4) it compiles ...

…lephase still needs testing

…t compositionalMultiphaseFlow/soreideWhitson/1D_100cells/1D_benchmark.xml

… to fixed and inputFiles regenerated

…d convert compositionalMultiphaseFlow/soreideWhitson/1D_100cells/1D_benchmark.xml" This reverts commit fa61ab1.

…well initialization vagaries assoiated with useSurfaceConditions=0

…for resvol constraint, compositional only

…te well state

joshua-white · 2026-07-15T00:20:57Z

      // set the reference well element where the BHP control is applied
-      wellControls.setReferenceGravityCoef( refElev * gravVector[2] );
+      wellControls.setReferenceGravityCoef( refElev * gravVector[2] );  // tjb remove


actually remove?

victorapm · 2026-07-16T20:20:37Z

+                                                                 rankNegPressureIds,
+                                                                 rankNegDensityIds,
+                                                                 rankTotalNegDensityIds );
+      globalCheck = MpiWrapper::min( localCheck );


Should we move this MPI call outside the for loop (for performance reasons)?

victorapm · 2026-07-16T20:24:11Z

+        for( integer i=0; i<numNorm; i++ )
+        {
+          localResidualNorm[i] = 0.0;
+        }


If I'm reading it correctly, there might be an issue when an open unconverged well is followed by a closed well. This branch resets the manager's accumulated maximum residual to zero, while I think closed wells should contribute nothing without clearing residuals already calculated for other wells.

I think optionizing the wellResidualNorm on open wells and then max the localResidualNorm sounds safer

if( wellControls.isWellOpen() ) { wellResidualNorm = wellControls.calculateLocalWellResidualNorm( time_n, dt, m_nonlinearSolverParameters, subRegion, dofManager, localRhs ); for( integer i = 0; i < numNorm; i++ ) { localResidualNorm[i] = LvArray::math::max( localResidualNorm[i], wellResidualNorm[i] ); } }

victorapm · 2026-07-16T21:33:01Z

+ENUM_STRINGS( WellControls_Control,
              "BHP",
              "phaseVolRate",
              "totalVolRate",
              "massRate",
              "uninitialized" );


The enum ConstraintTypeId has 7 values, while we are listing only five here. Should we update this list?

victorapm · 2026-07-16T23:14:30Z

+      // Assumes useMass is true
+      currentMassRate = currentTotalRate;


Should we add a check here on useMass instead of assuming it's true?

victorapm · 2026-07-16T23:16:02Z

+}
+void WellManager::initializePostSubGroups()
+{
+#if 0


Just curious why this code block is ifdef'ed out

victorapm · 2026-07-17T00:31:02Z

+  // 13.4) Make sure there is at least one non-BHP constraint
+  bool const rate_match_found = std::any_of( constraints.begin(), constraints.end(), [&bhp_type]( const auto & constraint_tuple )
+  {
+    return std::get< 1 >( constraint_tuple ) != bhp_type;
+  } );
+  GEOS_THROW_IF( !rate_match_found,
+                 GEOS_FMT( "Missing rate constraint for {} well {}",
+                           (isProducerWell ? "producer" : "injector"), getName() ),


if a well has phaseVolRate while defining only a total-volume rate constraint, this check would flag things are good, but i guess initialization would hit issues later.

victorapm · 2026-07-17T00:34:40Z

+  string const bhp_type = isProducerWell ? MinimumBHPConstraint::catalogName() : MaximumBHPConstraint::catalogName();
+  bool const no_match_found = std::none_of( constraints.begin(), constraints.end(), [&bhp_type]( const auto & constraint_tuple )
+  {
+    return std::get< 1 >( constraint_tuple ) == bhp_type;
+  } );
+  GEOS_THROW_IF( no_match_found,


When the sole BHP or rate constraint has constraintActive=0, these existence checks still accept it, while getBHPConstraint and the rate-list helpers filter it out. Subsequent selection either dereferences null or erases an end() iterator because the current constraint is absent from the active list.

victorapm · 2026-07-17T01:08:25Z

+    if( !hasZeroRate )
+    {
+      for( auto const * constraint : getRateConstraints() )
+      {
+        if( isZero( constraint->getConstraintValue( currentTime ) ) )
+        {
+          hasZeroRate = true;
+        }
+      }
+    }


When one active secondary rate schedule evaluates to zero while another rate target remains positive, this loop closes the entire well on the first zero (code below). Should we close the well only when every active rate target is zero?

victorapm · 2026-07-17T01:12:45Z

+      well.setThermal( isThermal() );
+      well.setUseMass( m_useMass );
+      well.registerWellDataOnMesh( subRegion );


Shouldn't this also be setting useTotalMassEquation and allowLocalCompDensityChopping?

victorapm · 2026-07-17T14:19:24Z

+  forSubGroups< InjectionConstraint< PhaseVolumeRateConstraint >, ProductionConstraint< PhaseVolumeRateConstraint > >( [&]( auto & constraint )
+  {
+    constraint.validatePhaseType( fluid );
+  } );


When a liquid-rate constraint becomes limiting, only phase-volume constraints have their phase names validated here. The liquid constraint's m_phaseIndices remains empty, which I believe could lead to issues during assembly (zero Jacobian).

victorapm · 2026-07-17T14:23:48Z

+{
+  real64 const currentValue = currentConstraint.liquidRate();
+  real64 const constraintValue = this->getConstraintValue( currentTime );
+  return ( LvArray::math::abs( currentValue ) <= LvArray::math::abs( constraintValue ) );


If I read correctly, for a production liquid-rate maximum, this comparison marks values below the target as violations and accepts values above it. This is the reverse of the other maximum-rate constraints and switches control incorrectly.

jafranc · 2026-07-17T16:02:37Z

+   * @param domain the physical domain object
+   */
+  void chopNegativeDensities( DomainPartition & domain );
+#if 0


should it be gone ?

jafranc · 2026-07-17T16:12:11Z

+  m_useMass( false ),
+  m_useTotalMassEquation( 1 ),
+  m_isThermal( 0 ),
+  m_isCompositional( true ),


should'nt it be an integer and registered for restart ? what happen if we want to restard a non-compositional well ?

jafranc · 2026-07-17T16:26:25Z

+  array1d< real64 > localResidualNorm, wellResidalNorm;
+  array1d< real64 > localResidualNormalizer;
+
+  if( isThermal() )


Sticking to conventionm shouldn't it be set by NDOF or so ?

jafranc · 2026-07-17T16:32:38Z

+          if( wellResidalNorm[i] > localResidualNorm[i] )
+          {
+            localResidualNorm[i] = wellResidalNorm[i];
+          }


Suggested change

if( wellResidalNorm[i] > localResidualNorm[i] )

{

localResidualNorm[i] = wellResidalNorm[i];

}

localResidualNorm[i] = LvArray::math::max( localResidualNorm[i], wellResidualNorm[i] );

jafranc · 2026-07-17T16:33:56Z

+      // step 1: compute the norm in the subRegion
+      if( wellControls.isWellOpen( ) )
+      {
+        wellResidalNorm = wellControls.calculateLocalWellResidualNorm( time_n,


Suggested change

wellResidalNorm = wellControls.calculateLocalWellResidualNorm( time_n,

wellResidualNorm = wellControls.calculateLocalWellResidualNorm( time_n,

jafranc · 2026-07-17T16:46:40Z

+        for( integer i=0; i<numNorm; i++ )
+        {
+          localResidualNorm[i] = 0.0;
+        }


I think optionizing the wellResidualNorm on open wells and then max the localResidualNorm sounds safer

if( wellControls.isWellOpen() ) { wellResidualNorm = wellControls.calculateLocalWellResidualNorm( time_n, dt, m_nonlinearSolverParameters, subRegion, dofManager, localRhs ); for( integer i = 0; i < numNorm; i++ ) { localResidualNorm[i] = LvArray::math::max( localResidualNorm[i], wellResidualNorm[i] ); } }

tjb-ltk and others added 30 commits September 29, 2025 16:50

1) new constraint schema, 2) resturcture well jac assem loops

f9db9fd

fix merge errors

a9cd276

fix integrate errors

c59b1a6

rel mode fix

48aabd4

1) remove old code, 2) test different constraint org, 3) singlephase …

a85ee1d

…broken

1) SinglePhaseWell adopts new constraint convention, 2) Change constr…

8366922

…aint class hierachy , 3) Remove all old constraint handling methods/variables, 4) it compiles ...

fix constraint sign and normalizers

4f49bf1

1) compositional_multiphase_wells_1d.xml matches geos branch, 2) sing…

2e82b49

…lephase still needs testing

compressible_single_phase_wells_1d.xml matches geos branch

6281018

improve naming conventions for well constriant files and classes

205889f

fix reference to 'internal' is ambiguous error

b6ae5ef

1) well constraint schema script, 2) converted xmls, 3) Failed conver…

fa61ab1

…t compositionalMultiphaseFlow/soreideWhitson/1D_100cells/1D_benchmark.xml

schema

292c90f

1)fixes for catalogue registration, 2) xml convert script still needs…

d4eb504

… to fixed and inputFiles regenerated

fix liq rate schema

bd3d0cf

Revert "1) well constraint schema script, 2) converted xmls, 3) Faile…

723c2bc

…d convert compositionalMultiphaseFlow/soreideWhitson/1D_100cells/1D_benchmark.xml" This reverts commit fa61ab1.

xml with new constraint schema, try 2

e39ccf1

merge devel

d0fcf79

remove extra ;

8b21d69

remove :

470a7ff

fix constraint schema, phase name on totalvol constraint not allowed

b8a90f4

update well constraint schema

94fc9ad

1) expose method for unit tests, 2) update unittests accordingly

17047bd

schema change

c6a349f

1) remove unused vars , add missing constriant type in selector

1c03220

reintroduce control keyword option to be consistent with previous in …

b5f6c5d

…well initialization vagaries assoiated with useSurfaceConditions=0

updated schema

c43657e

add back control keyword

e78ef18

add control schema to unit tests

b000ed4

uncrustify should be automated

a319e0e

tjb-ltk added 6 commits June 22, 2026 08:12

Add comments to show activation of various derivative tests

a83f585

remove duplicate entry

a43b83e

Merge branch 'develop' into feature/byer3/wm_we

8fde458

vscode merge isnt always crustify compliant

37c429b

missing file for well schema converter

08f17b3

minor cleanup

19f2305

joshua-white reviewed Jul 9, 2026

View reviewed changes

Comment thread src/coreComponents/integrationTests/wellsTests/testThermalInjWell.cpp

joshua-white reviewed Jul 9, 2026

View reviewed changes

Comment thread ...coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp Outdated

joshua-white reviewed Jul 9, 2026

View reviewed changes

Comment thread src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/SinglePhaseWellKernels.cpp Outdated

tjb-ltk added 5 commits July 10, 2026 10:34

Merge branch 'develop' into feature/byer3/wm_we

5e2935e

more devel merge induced changes

ade7016

1) pr review changes, 2) code refactor to support regionstats needed …

005081f

…for resvol constraint, compositional only

region statistics interface fixes 1) cut timestep , 2) estimator upda…

e0f350a

…te well state

schema update

c0390d6

joshua-white reviewed Jul 15, 2026

View reviewed changes

joshua-white approved these changes Jul 15, 2026

View reviewed changes

integrated test update due to schema and data structure changes

ea74e49

herve-gross approved these changes Jul 16, 2026

View reviewed changes

victorapm reviewed Jul 16, 2026

View reviewed changes

victorapm reviewed Jul 17, 2026

View reviewed changes

jafranc reviewed Jul 17, 2026

View reviewed changes

		// Assumes useMass is true
		currentMassRate = currentTotalRate;

	wellResidalNorm = wellControls.calculateLocalWellResidualNorm( time_n,
	wellResidualNorm = wellControls.calculateLocalWellResidualNorm( time_n,

Uh oh!

Conversation

tjb-ltk commented Feb 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Overview

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

8 participants

tjb-ltk commented Feb 13, 2026 •

edited

Loading