feat: well time step selector based on rates/bhp tables and clarify well rates logic

src/coreComponents/functions/TableFunction.cpp

@@ -196,6 +196,11 @@ 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
+{
+  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 ...
+     * @param[in] input vector of input value
+     * @param[in] dim table dimension
+     * @param[in] interpolationMethod interpolation method
+     * @return coordinate value
+     */
+    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

@@ -302,7 +302,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)
@@ -332,7 +332,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 )
@@ -363,15 +363,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 nextDtNewton = setNextDtBasedOnIterNumber( currentDt );
+  if( m_nonlinearSolverParameters.getLogLevel() > 0 && nextDtNewton < LvArray::NumericLimits< real64 >::max )
+    GEOS_LOG_RANK_0( GEOS_FMT( "{}: next time step based on number of iterations = {}", getName(), nextDtNewton ));
   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 ) )
@@ -428,7 +429,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();
@@ -458,15 +459,6 @@ real64 PhysicsSolverBase::setNextDtBasedOnNewtonIter( real64 const & currentDt )
   return nextDt;
 }
 
-
-real64 PhysicsSolverBase::setNextDtBasedOnCFL( const geos::real64 & currentDt, geos::DomainPartition & domain )
-{
-  GEOS_UNUSED_VAR( currentDt, domain );
-  return LvArray::NumericLimits< real64 >::max;       // i.e., not implemented
-}
-
-
-
 real64 PhysicsSolverBase::linearImplicitStep( real64 const & time_n,
                                               real64 const & dt,
                                               integer const GEOS_UNUSED_PARAM( cycleNumber ),

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
@@ -239,17 +241,6 @@ class PhysicsSolverBase : public ExecutableGroup
   virtual real64 setNextDtBasedOnStateChange( real64 const & currentDt,
                                               DomainPartition & domain );
 
-  /**
-   * @brief function to set the next dt based on state change
-   * @param [in]  currentDt the current time step size
-   * @param[in] domain the domain object
-   * @return the prescribed time step size
-   */
-  virtual real64 setNextDtBasedOnCFL( real64 const & currentDt,
-                                      DomainPartition & domain );
-
-
-
   /**
    * @brief Entry function for an explicit time integration step
    * @param time_n time at the beginning of the step

src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp

@@ -2355,13 +2355,6 @@ bool SolidMechanicsLagrangeContact::isFractureAllInStickCondition( DomainPartiti
   return ( ( numNewSlip + numSlip + numOpen ) == 0 );
 }
 
-real64 SolidMechanicsLagrangeContact::setNextDt( real64 const & currentDt,
-                                                 DomainPartition & domain )
-{
-  GEOS_UNUSED_VAR( domain );
-  return currentDt;
-}
-
 REGISTER_CATALOG_ENTRY( PhysicsSolverBase, SolidMechanicsLagrangeContact, string const &, Group * const )
 
 } /* namespace geos */

src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.hpp

@@ -101,10 +101,6 @@ class SolidMechanicsLagrangeContact : public ContactSolverBase
   virtual void
   resetStateToBeginningOfStep( DomainPartition & domain ) override;
 
-  virtual real64
-  setNextDt( real64 const & currentDt,
-             DomainPartition & domain ) override;
-
   void updateState( DomainPartition & domain ) override final;
 
   void assembleContact( DomainPartition & domain,

src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp

@@ -2653,11 +2653,12 @@ bool CompositionalMultiphaseBase::checkSequentialSolutionIncrements( DomainParti
   return isConverged && (m_sequentialCompDensChange < m_maxSequentialCompDensChange);
 }
 
-real64 CompositionalMultiphaseBase::setNextDt( const geos::real64 & currentDt, geos::DomainPartition & domain )
+real64 CompositionalMultiphaseBase::setNextDt( real64 const & currentTime,
+                                               const real64 & currentDt,
+                                               DomainPartition & domain )
 {
-
   if( m_targetFlowCFL < 0 )
-    return PhysicsSolverBase::setNextDt( currentDt, domain );
+    return PhysicsSolverBase::setNextDt( currentTime, currentDt, domain );
   else
     return setNextDtBasedOnCFL( currentDt, domain );
 }

src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp

@@ -362,22 +362,24 @@ class CompositionalMultiphaseBase : public FlowSolverBase
    */
   void chopNegativeDensities( DomainPartition & domain );
 
-  virtual real64 setNextDtBasedOnStateChange( real64 const & currentDt,
-                                              DomainPartition & domain ) override;
-
   void computeCFLNumbers( DomainPartition & domain, real64 const & dt, real64 & maxPhaseCFL, real64 & maxCompCFL );
 
   /**
    * @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,
-                    DomainPartition & domain ) override;
+  virtual real64 setNextDt( real64 const & currentTime,
+                            real64 const & currentDt,
+                            DomainPartition & domain ) override;
+
+  virtual real64 setNextDtBasedOnStateChange( real64 const & currentDt,
+                                              DomainPartition & domain ) override;
 
-  virtual real64 setNextDtBasedOnCFL( real64 const & currentDt,
-                                      DomainPartition & domain ) override;
+  real64 setNextDtBasedOnCFL( real64 const & currentDt,
+                              DomainPartition & domain );
 
   virtual void initializePostInitialConditionsPreSubGroups() override;