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 );
 }
 
+std::vector< real64 > TableFunction::getCoord( real64 const * const input, InterpolationType interpolationMethod ) const
+{
+  return m_kernelWrapper.getCoord( input, 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,16 @@ class TableFunction : public FunctionBase
     real64
     interpolateRound( IN_ARRAY const & input ) const;
 
+    /**
+     * @brief ...
+     * @param[in] input vector of input value
+     * @return coordinate value
+     */
+    template< typename IN_ARRAY >
+    GEOS_HOST_DEVICE
+    std::vector< real64 >
+    getCoord( IN_ARRAY const & input, InterpolationType interpolationMethod ) const;
+
     /**
      * @brief Interpolate in the table with derivatives using linear method.
      * @param[in] input vector of input value
@@ -240,6 +250,14 @@ 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 interpolationMethod the interpolation method
+   * @return the coordinates
+   */
+  std::vector< real64 > getCoord( real64 const * const input, 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 +572,64 @@ TableFunction::KernelWrapper::interpolateRound( IN_ARRAY const & input ) const
   return m_values[tableIndex];
 }
 
+template< typename IN_ARRAY >
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+std::vector< real64 >
+TableFunction::KernelWrapper::getCoord( IN_ARRAY const & input, InterpolationType interpolationMethod ) const
+{
+  integer const numDimensions = LvArray::integerConversion< integer >( m_coordinates.size() );
+  std::vector< real64 > ret( numDimensions, 0.0 );
+
+  // Determine the index to the nearest table entry
+  localIndex subIndex;
+  for( integer dim = 0; dim < numDimensions; ++dim )
+  {
+    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
+    ret[dim] = coords[subIndex];
+  }
+
+  return ret;
+}
+
 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,50 +363,51 @@ void PhysicsSolverBase::logEndOfCycleInformation( integer const cycleNumber,
   GEOS_LOG_RANK_0( "------------------------------------------------------------------\n" );
 }
 
-real64 PhysicsSolverBase::setNextDt( real64 const & currentDt,
+real64 PhysicsSolverBase::setNextDt( real64 const & GEOS_UNUSED_PARAM( currentTime ),
+                                     real64 const & lastDt,
                                      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 nextDtStateChange = setNextDtBasedOnStateChange( currentDt, domain );
-  if( m_nonlinearSolverParameters.getLogLevel() > 0 )
+  real64 const nextDtNewton = setNextDtBasedOnIterNumber( lastDt );
+  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( lastDt, domain );
+  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 ) )
   {
-    GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: time-step size will be kept the same since it's been {} cycles since last cut.",
+    GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: time-step size should be kept the same since it's been {} cycles since last cut",
                                                              getName(), m_numTimestepsSinceLastDtCut ) );
-    return currentDt;
+    return lastDt;
   }
 
   if( nextDtNewton < nextDtStateChange )      // time step size decided based on convergence
   {
-    if( nextDtNewton > currentDt )
+    if( nextDtNewton > lastDt )
     {
-      GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: time-step required will be increased based on number of iterations.",
+      GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: time-step required should be increased based on number of iterations",
                                                                getName() ) );
     }
-    else if( nextDtNewton < currentDt )
+    else if( nextDtNewton < lastDt )
     {
-      GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: time-step required will be decreased based on number of iterations.",
+      GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: time-step required should be decreased based on number of iterations",
                                                                getName() ) );
     }
     else
     {
-      GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: time-step required will be kept the same based on number of iterations.",
+      GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: time-step required should be kept the same based on number of iterations",
                                                                getName() ) );
     }
   }
   else         // time step size decided based on state change
   {
-    if( nextDtStateChange > currentDt )
+    if( nextDtStateChange > lastDt )
     {
       GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: time-step required will be increased based on state change.",
                                                                getName()));
     }
-    else if( nextDtStateChange < currentDt )
+    else if( nextDtStateChange < lastDt )
     {
       GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::TimeStep, GEOS_FMT( "{}: time-step required will be decreased based on state change.",
                                                                getName()));
@@ -421,14 +422,14 @@ real64 PhysicsSolverBase::setNextDt( real64 const & currentDt,
   return std::min( nextDtNewton, nextDtStateChange );
 }
 
-real64 PhysicsSolverBase::setNextDtBasedOnStateChange( real64 const & currentDt,
+real64 PhysicsSolverBase::setNextDtBasedOnStateChange( real64 const & lastDt,
                                                        DomainPartition & domain )
 {
-  GEOS_UNUSED_VAR( currentDt, domain );
+  GEOS_UNUSED_VAR( lastDt, domain );
   return LvArray::NumericLimits< real64 >::max; // i.e., not implemented
 }
 
-real64 PhysicsSolverBase::setNextDtBasedOnNewtonIter( real64 const & currentDt )
+real64 PhysicsSolverBase::setNextDtBasedOnIterNumber( real64 const & lastDt )
 {
   integer & newtonIter = m_nonlinearSolverParameters.m_numNewtonIterations;
   integer const iterDecreaseLimit = m_nonlinearSolverParameters.timeStepDecreaseIterLimit();
@@ -438,35 +439,26 @@ real64 PhysicsSolverBase::setNextDtBasedOnNewtonIter( real64 const & currentDt )
   if( newtonIter < iterIncreaseLimit )
   {
     // Easy convergence, let's increase the time-step.
-    nextDt = currentDt * m_nonlinearSolverParameters.timeStepIncreaseFactor();
+    nextDt = lastDt * m_nonlinearSolverParameters.timeStepIncreaseFactor();
     if( m_nonlinearSolverParameters.getLogLevel() > 0 )
       GEOS_LOG_RANK_0( GEOS_FMT( "{}: number of iterations = {} is less than {}, next time step = {} (increase)", getName(), newtonIter, iterIncreaseLimit, nextDt ));
   }
   else if( newtonIter > iterDecreaseLimit )
   {
     // Tough convergence let us make the time-step smaller!
-    nextDt = currentDt * m_nonlinearSolverParameters.timeStepDecreaseFactor();
+    nextDt = lastDt * m_nonlinearSolverParameters.timeStepDecreaseFactor();
     if( m_nonlinearSolverParameters.getLogLevel() > 0 )
       GEOS_LOG_RANK_0( GEOS_FMT( "{}: number of iterations = {} is more than {}, next time step = {} (decrease)", getName(), newtonIter, iterDecreaseLimit, nextDt ));
   }
   else
   {
-    nextDt = currentDt;
+    nextDt = lastDt;
     if( m_nonlinearSolverParameters.getLogLevel() > 0 )
       GEOS_LOG_RANK_0( GEOS_FMT( "{}: number of iterations = {} is between {} and {}, next time step = {} (no change)", getName(), newtonIter, iterIncreaseLimit, iterDecreaseLimit, nextDt ));
   }
   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,40 +216,31 @@ class PhysicsSolverBase : public ExecutableGroup
 
   /**
    * @brief function to set the next time step size
-   * @param[in] currentDt the current time step size
+   * @param[in] time the current time
+   * @param[in] lastDt the last 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 & time,
+                            real64 const & lastDt,
                             DomainPartition & domain );
 
   /**
-   * @brief function to set the next time step size based on Newton convergence
-   * @param[in] currentDt the current time step size
+   * @brief function to set the next time step size based on convergence
+   * @param[in] lastDt the last time step size
    * @return the prescribed time step size
    */
-  virtual real64 setNextDtBasedOnNewtonIter( real64 const & currentDt );
+  virtual real64 setNextDtBasedOnIterNumber( real64 const & lastDt );
 
   /**
    * @brief function to set the next dt based on state change
-   * @param [in]  currentDt the current time step size
+   * @param [in]  lastDt the last time step size
    * @param[in] domain the domain object
    * @return the prescribed time step size
    */
-  virtual real64 setNextDtBasedOnStateChange( real64 const & currentDt,
+  virtual real64 setNextDtBasedOnStateChange( real64 const & lastDt,
                                               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

@@ -2356,13 +2356,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,