feat: well time step selector based on rates/bhp tables and clarify well rates logic (#3427)

Author: paveltomin

.integrated_tests.yaml

@@ -1,6 +1,6 @@
 baselines:
   bucket: geosx
-  baseline: integratedTests/baseline_integratedTests-pr3485-10612-832f30e
+  baseline: integratedTests/baseline_integratedTests-pr3427-10674-34d1de8
 
 allow_fail:
   all: ''

BASELINE_NOTES.md

@@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines.
 Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining.
 These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD).
 
+PR #3427 (2024-03-10)
+=====================
+Well time step selector based on rates/bhp tables and clarify well rates logic.
+
 PR #3485 (2024-03-09)
 =====================
 Use mass and energy consistently for single phase solvers.

src/coreComponents/codingUtilities/Utilities.hpp

@@ -59,6 +59,20 @@ bool isZero( T const val, T const tol = LvArray::NumericLimits< T >::epsilon )
   return -tol <= val && val <= tol;
 }
 
+template< typename ARRAY_TYPE >
+GEOS_FORCE_INLINE GEOS_HOST_DEVICE
+bool hasNonZero( ARRAY_TYPE const & array )
+{
+  for( auto it = array.begin(); it != array.end() ; ++it )
+  {
+    if( !isZero( *it ) )
+    {
+      return true;
+    }
+  }
+  return false;
+}
+
 template< typename T >
 GEOS_FORCE_INLINE GEOS_HOST_DEVICE constexpr
 bool isOdd( T x )

src/coreComponents/functions/TableFunction.cpp

@@ -196,6 +196,15 @@ real64 TableFunction::evaluate( real64 const * const input ) const
   return m_kernelWrapper.compute( input );
 }
 
+real64 TableFunction::getCoord( real64 const * const input, localIndex const dim, InterpolationType interpolationMethod ) const
+{
+  GEOS_ASSERT_MSG( interpolationMethod != InterpolationType::Linear,
+                   GEOS_FMT( "{}: TableFunction::getCoord should not be called with {} interpolation method",
+                             getDataContext(), EnumStrings< InterpolationType >::toString( interpolationMethod )));
+  GEOS_ASSERT( dim >= 0 && dim < m_coordinates.size() );
+  return m_kernelWrapper.getCoord( input, dim, interpolationMethod );
+}
+
 TableFunction::KernelWrapper::KernelWrapper( InterpolationType const interpolationMethod,
                                              ArrayOfArraysView< real64 const > const & coordinates,
                                              arrayView1d< real64 const > const & values )

src/coreComponents/functions/TableFunction.hpp

@@ -173,6 +173,18 @@ class TableFunction : public FunctionBase
     real64
     interpolateRound( IN_ARRAY const & input ) const;
 
+    /**
+     * @brief Method to get coordinates
+     * @param input a scalar input
+     * @param dim the table dimension
+     * @param interpolationMethod the interpolation method
+     * @return the coordinate
+     */
+    template< typename IN_ARRAY >
+    GEOS_HOST_DEVICE
+    real64
+    getCoord( IN_ARRAY const & input, localIndex dim, InterpolationType interpolationMethod ) const;
+
     /**
      * @brief Interpolate in the table with derivatives using linear method.
      * @param[in] input vector of input value
@@ -240,6 +252,15 @@ class TableFunction : public FunctionBase
    */
   virtual real64 evaluate( real64 const * const input ) const override final;
 
+  /**
+   * @brief Method to get coordinates
+   * @param input a scalar input
+   * @param dim the table dimension
+   * @param interpolationMethod the interpolation method
+   * @return the coordinate
+   */
+  real64 getCoord( real64 const * const input, localIndex dim, InterpolationType interpolationMethod ) const;
+
   /**
    * @brief Check if the given coordinate is in the bounds of the table coordinates in the
    * specified dimension, throw an exception otherwise.
@@ -554,6 +575,57 @@ TableFunction::KernelWrapper::interpolateRound( IN_ARRAY const & input ) const
   return m_values[tableIndex];
 }
 
+template< typename IN_ARRAY >
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+real64
+TableFunction::KernelWrapper::getCoord( IN_ARRAY const & input, localIndex const dim, InterpolationType interpolationMethod ) const
+{
+  // Determine the index to the nearest table entry
+  localIndex subIndex;
+  arraySlice1d< real64 const > const coords = m_coordinates[dim];
+  // Determine the index along each table axis
+  if( input[dim] <= coords[0] )
+  {
+    // Coordinate is to the left of the table axis
+    subIndex = 0;
+  }
+  else if( input[dim] >= coords[coords.size() - 1] )
+  {
+    // Coordinate is to the right of the table axis
+    subIndex = coords.size() - 1;
+  }
+  else
+  {
+    // Coordinate is within the table axis
+    // Note: find() will return the index of the upper table vertex
+    auto const lower = LvArray::sortedArrayManipulation::find( coords.begin(), coords.size(), input[dim] );
+    subIndex = LvArray::integerConversion< localIndex >( lower );
+
+    // Interpolation types:
+    //   - Nearest returns the value of the closest table vertex
+    //   - Upper returns the value of the next table vertex
+    //   - Lower returns the value of the previous table vertex
+    if( interpolationMethod == TableFunction::InterpolationType::Nearest )
+    {
+      if( ( input[dim] - coords[subIndex - 1]) <= ( coords[subIndex] - input[dim]) )
+      {
+        --subIndex;
+      }
+    }
+    else if( interpolationMethod == TableFunction::InterpolationType::Lower )
+    {
+      if( subIndex > 0 )
+      {
+        --subIndex;
+      }
+    }
+  }
+
+  // Retrieve the nearest coordinate
+  return coords[subIndex];
+}
+
 template< typename IN_ARRAY, typename OUT_ARRAY >
 GEOS_HOST_DEVICE
 GEOS_FORCE_INLINE

src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp

@@ -303,7 +303,7 @@ bool PhysicsSolverBase::execute( real64 const time_n,
 
     if( dtRemaining > 0.0 )
     {
-      nextDt = setNextDt( dtAccepted, domain );
+      nextDt = setNextDt( time_n + (dt - dtRemaining), dtAccepted, domain );
       if( nextDt < dtRemaining )
       {
         // better to do two equal steps than one big and one small (even potentially tiny)
@@ -333,7 +333,7 @@ bool PhysicsSolverBase::execute( real64 const time_n,
                                                         " reached. Consider increasing maxSubSteps." );
 
   // Decide what to do with the next Dt for the event running the solver.
-  m_nextDt = setNextDt( nextDt, domain );
+  m_nextDt = setNextDt( time_n + dt, nextDt, domain );
 
   // Increase counter to indicate how many cycles since the last timestep cut
   if( m_numTimestepsSinceLastDtCut >= 0 )
@@ -364,15 +364,16 @@ void PhysicsSolverBase::logEndOfCycleInformation( integer const cycleNumber,
   GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, "------------------------------------------------------------------\n" );
 }
 
-real64 PhysicsSolverBase::setNextDt( real64 const & currentDt,
+real64 PhysicsSolverBase::setNextDt( real64 const & GEOS_UNUSED_PARAM( currentTime ),
+                                     real64 const & currentDt,
                                      DomainPartition & domain )
 {
   integer const minTimeStepIncreaseInterval = m_nonlinearSolverParameters.minTimeStepIncreaseInterval();
-  real64 const nextDtNewton = setNextDtBasedOnNewtonIter( currentDt );
-  if( m_nonlinearSolverParameters.getLogLevel() > 0 )
-    GEOS_LOG_RANK_0( GEOS_FMT( "{}: next time step based on Newton iterations = {}", getName(), nextDtNewton ));
+  real64 const nextDtIter = setNextDtBasedOnIterNumber( currentDt );
+  if( m_nonlinearSolverParameters.getLogLevel() > 0 && nextDtIter < LvArray::NumericLimits< real64 >::max )
+    GEOS_LOG_RANK_0( GEOS_FMT( "{}: next time step based on number of iterations = {}", getName(), nextDtIter ));
   real64 const nextDtStateChange = setNextDtBasedOnStateChange( currentDt, domain );
-  if( m_nonlinearSolverParameters.getLogLevel() > 0 )
+  if( m_nonlinearSolverParameters.getLogLevel() > 0 && nextDtStateChange < LvArray::NumericLimits< real64 >::max )
     GEOS_LOG_RANK_0( GEOS_FMT( "{}: next time step based on state change = {}", getName(), nextDtStateChange ));
 
   if( ( m_numTimestepsSinceLastDtCut >= 0 ) && ( m_numTimestepsSinceLastDtCut < minTimeStepIncreaseInterval ) )
@@ -382,14 +383,14 @@ real64 PhysicsSolverBase::setNextDt( real64 const & currentDt,
     return currentDt;
   }
 
-  if( nextDtNewton < nextDtStateChange )      // time step size decided based on convergence
+  if( nextDtIter < nextDtStateChange )      // time step size decided based on convergence
   {
-    if( nextDtNewton > currentDt )
+    if( nextDtIter > currentDt )
     {
       GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: time-step required will be increased based on number of iterations.",
                                                                getName() ) );
     }
-    else if( nextDtNewton < currentDt )
+    else if( nextDtIter < currentDt )
     {
       GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: time-step required will be decreased based on number of iterations.",
                                                                getName() ) );
@@ -419,7 +420,7 @@ real64 PhysicsSolverBase::setNextDt( real64 const & currentDt,
     }
   }
 
-  return std::min( nextDtNewton, nextDtStateChange );
+  return std::min( nextDtIter, nextDtStateChange );
 }
 
 real64 PhysicsSolverBase::setNextDtBasedOnStateChange( real64 const & currentDt,
@@ -429,7 +430,7 @@ real64 PhysicsSolverBase::setNextDtBasedOnStateChange( real64 const & currentDt,
   return LvArray::NumericLimits< real64 >::max; // i.e., not implemented
 }
 
-real64 PhysicsSolverBase::setNextDtBasedOnNewtonIter( real64 const & currentDt )
+real64 PhysicsSolverBase::setNextDtBasedOnIterNumber( real64 const & currentDt )
 {
   integer & newtonIter = m_nonlinearSolverParameters.m_numNewtonIterations;
   integer const iterDecreaseLimit = m_nonlinearSolverParameters.timeStepDecreaseIterLimit();

src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp

@@ -216,19 +216,21 @@ class PhysicsSolverBase : public ExecutableGroup
 
   /**
    * @brief function to set the next time step size
+   * @param[in] currentTime the current time
    * @param[in] currentDt the current time step size
    * @param[in] domain the domain object
    * @return the prescribed time step size
    */
-  virtual real64 setNextDt( real64 const & currentDt,
+  virtual real64 setNextDt( real64 const & currentTime,
+                            real64 const & currentDt,
                             DomainPartition & domain );
 
   /**
-   * @brief function to set the next time step size based on Newton convergence
+   * @brief function to set the next time step size based on convergence
    * @param[in] currentDt the current time step size
    * @return the prescribed time step size
    */
-  virtual real64 setNextDtBasedOnNewtonIter( real64 const & currentDt );
+  virtual real64 setNextDtBasedOnIterNumber( real64 const & currentDt );
 
   /**
    * @brief function to set the next dt based on state change

src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp

@@ -1554,10 +1554,12 @@ void CompositionalMultiphaseFVM::assembleHydrofracFluxTerms( real64 const GEOS_U
   } );
 }
 
-real64 CompositionalMultiphaseFVM::setNextDt( const geos::real64 & currentDt, geos::DomainPartition & domain )
+real64 CompositionalMultiphaseFVM::setNextDt( real64 const & currentTime,
+                                              real64 const & currentDt,
+                                              DomainPartition & domain )
 {
   if( m_targetFlowCFL < 0 )
-    return CompositionalMultiphaseBase::setNextDt( currentDt, domain );
+    return CompositionalMultiphaseBase::setNextDt( currentTime, currentDt, domain );
   else
     return setNextDtBasedOnCFL( currentDt, domain );
 }

src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp

@@ -195,11 +195,13 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase
 
   /**
    * @brief function to set the next time step size
+   * @param[in] currentTime the current time
    * @param[in] currentDt the current time step size
    * @param[in] domain the domain object
    * @return the prescribed time step size
    */
-  real64 setNextDt( real64 const & currentDt,
+  real64 setNextDt( real64 const & currentTime,
+                    real64 const & currentDt,
                     DomainPartition & domain ) override;
 
   struct viewKeyStruct : CompositionalMultiphaseBase::viewKeyStruct