Fix performance regression of mesh map generation (#1847)

src/coreComponents/codingUtilities/Utilities.hpp

@@ -118,6 +118,69 @@ T_VALUE softMapLookup( mapBase< T_KEY, T_VALUE, SORTED > const & theMap,
   return rvalue;
 }
 
+/**
+ * @brief Call a user function on sub-ranges of repeating values.
+ * @tparam ITER type of range iterator
+ * @tparam FUNC type of user function
+ * @tparam COMP type of comparison function
+ * @param first iterator to start of the input range
+ * @param last iterator past the end of the input range
+ * @param func the function to call, must be callable with a pair of iterators
+ *
+ * User function will be called once per consecutive group of equal values, as defined
+ * by @p comp, with a pair of iterators to the beginning and past the end of each group.
+ * For example, given an input [1,1,2,2,2,3,3,1], @p func will be called with iterators
+ * to sub-ranges [1,1], [2,2,2], [3,3], [1].
+ *
+ * @note @p comp is an equality comparison, not a less-than type predicate used in sorting.
+ *       For a range sorted with some predicate P, one can use comp(x,y) = !(P(x,y) || P(y,x)),
+ *       but in most cases default value (std::equal_to<>) should be fine.
+ */
+template< typename ITER, typename FUNC, typename COMP = std::equal_to<> >
+void forEqualRanges( ITER first, ITER const last, FUNC && func, COMP && comp = {} )
+{
+  using R = typename std::iterator_traits< ITER >::reference;
+  while( first != last )
+  {
+    R curr = *first;
+    auto const pred = [&]( R v ) { return comp( curr, v ); };
+    ITER const next = std::find_if_not( std::next( first ), last, pred );
+    func( first, next );
+    first = next;
+  }
+}
+
+/**
+ * @brief Execute a user function on unique values in a range.
+ * @tparam ITER type of range iterator
+ * @tparam FUNC type of user function
+ * @tparam COMP type of comparison function
+ * @param first iterator to start of the range
+ * @param last iterator past the end of the range
+ * @param func the function to call, must be callable with value and size (as int)
+ *
+ * User function will be called once per consecutive group of equal values, as defined
+ * by @p comp, with a reference to one of the values from each group and the group size.
+ * If the range is (partially) ordered w.r.t. to a predicate compatible with @p comp,
+ * the function will be called once per unique value in the entire range.
+ * For example, given an input [a,a,b,b,b,c,c,a], @p func will be called with
+ * (a,2), (b,3), (c,2), (a,1), where a,b,c will be references to one of the
+ * corresponding elements in the input (which one exactly is unspecified).
+ *
+ * @note @p comp is an equality comparison, not a less-than type predicate used in sorting.
+ *       For a range sorted with some predicate P, one can use comp(x,y) = !(P(x,y) || P(y,x)),
+ *       but in most cases default value (std::equal_to<>) should be fine.
+ */
+template< typename ITER, typename FUNC, typename COMP = std::equal_to<> >
+void forUniqueValues( ITER first, ITER const last, FUNC && func, COMP && comp = {} )
+{
+  auto const f = [&func]( ITER r_first, ITER r_last )
+  {
+    func( *r_first, std::distance( r_first, r_last ) );
+  };
+  forEqualRanges( first, last, f, std::forward< COMP >( comp ) );
+}
+
 /**
  * @brief Perform lookup in a map of options and throw a user-friendly exception if not found.
  * @tparam KEY map key type

src/coreComponents/constitutive/fluid/PVTDriver.cpp

@@ -182,7 +182,7 @@ bool PVTDriver::execute( real64 const GEOSX_UNUSED_PARAM( time_n ),
 
 
 template< typename FLUID_TYPE >
-void PVTDriver::runTest( FLUID_TYPE & fluid, arrayView2d< real64 > & table )
+void PVTDriver::runTest( FLUID_TYPE & fluid, arrayView2d< real64 > const & table )
 {
   // get number of phases and components
 
@@ -221,17 +221,18 @@ void PVTDriver::runTest( FLUID_TYPE & fluid, arrayView2d< real64 > & table )
     compositionValues[0][i] /= sum;
   }
 
-  arraySlice1d< real64 const, compflow::USD_COMP - 1 > const composition = compositionValues[0];
+  arrayView2d< real64 const, compflow::USD_COMP > const composition = compositionValues;
 
   // perform fluid update using table (P,T) and save resulting total density, etc.
   // note: column indexing should be kept consistent with output file header below.
 
+  integer numSteps = m_numSteps;
   using ExecPolicy = typename FLUID_TYPE::exec_policy;
-  forAll< ExecPolicy >( 1, [=]  GEOSX_HOST_DEVICE ( integer const ei )
+  forAll< ExecPolicy >( 1, [=]  GEOSX_HOST_DEVICE ( localIndex const ei )
   {
-    for( integer n=0; n<=m_numSteps; ++n )
+    for( integer n = 0; n <= numSteps; ++n )
     {
-      kernelWrapper.update( ei, 0, table( n, PRES ), table( n, TEMP ), composition );
+      kernelWrapper.update( ei, 0, table( n, PRES ), table( n, TEMP ), composition[0] );
       table( n, TEMP+1 ) = totalDensity( ei, 0 );
 
       for( integer p=0; p<NP; ++p )

src/coreComponents/constitutive/fluid/PVTDriver.hpp

@@ -64,7 +64,7 @@ class PVTDriver : public TaskBase
    * @param table Table with input/output time history
    */
   template< typename FLUID_TYPE >
-  void runTest( FLUID_TYPE & fluid, arrayView2d< real64 > & table );
+  void runTest( FLUID_TYPE & fluid, arrayView2d< real64 > const & table );
 
   /**
    * @brief Ouput table to file for easy plotting

src/coreComponents/constitutive/solid/TriaxialDriver.cpp

@@ -157,17 +157,18 @@ void TriaxialDriver::postProcessInput()
 
 
 template< typename SOLID_TYPE >
-void TriaxialDriver::runStrainControlTest( SOLID_TYPE & solid, arrayView2d< real64 > & table )
+void TriaxialDriver::runStrainControlTest( SOLID_TYPE & solid, arrayView2d< real64 > const & table )
 {
   typename SOLID_TYPE::KernelWrapper updates = solid.createKernelUpdates();
+  localIndex const numSteps = m_numSteps;
 
   forAll< parallelDevicePolicy<> >( 1, [=]  GEOSX_HOST_DEVICE ( integer const ei )
   {
     real64 stress[6] = {};
     real64 strainIncrement[6] = {};
     real64 stiffness[6][6] = {{}};
 
-    for( integer n=1; n<=m_numSteps; ++n )
+    for( integer n = 1; n <= numSteps; ++n )
     {
       strainIncrement[0] = table( n, EPS0 )-table( n-1, EPS0 );
       strainIncrement[1] = table( n, EPS1 )-table( n-1, EPS1 );
@@ -187,9 +188,13 @@ void TriaxialDriver::runStrainControlTest( SOLID_TYPE & solid, arrayView2d< real
 
 
 template< typename SOLID_TYPE >
-void TriaxialDriver::runMixedControlTest( SOLID_TYPE & solid, arrayView2d< real64 > & table )
+void TriaxialDriver::runMixedControlTest( SOLID_TYPE & solid, arrayView2d< real64 > const & table )
 {
   typename SOLID_TYPE::KernelWrapper updates = solid.createKernelUpdates();
+  integer const numSteps = m_numSteps;
+  integer const maxIter = m_maxIter;
+  integer const maxCuts = m_maxCuts;
+  real64 const newtonTol = m_newtonTol;
 
   forAll< parallelDevicePolicy<> >( 1, [=]  GEOSX_HOST_DEVICE ( integer const ei )
   {
@@ -199,13 +204,13 @@ void TriaxialDriver::runMixedControlTest( SOLID_TYPE & solid, arrayView2d< real6
     real64 stiffness[6][6] = {{}};
 
     real64 scale = 0;
-    for( integer n=1; n<=m_numSteps; ++n )
+    for( integer n = 1; n <= numSteps; ++n )
     {
       scale += fabs( table( n, SIG0 )) + fabs( table( n, SIG1 )) + fabs( table( n, SIG2 ));
     }
-    scale = 3*m_numSteps / scale;
+    scale = 3 * numSteps / scale;
 
-    for( integer n=1; n<=m_numSteps; ++n )
+    for( integer n=1; n<=numSteps; ++n )
     {
       strainIncrement[0] = table( n, EPS0 )-table( n-1, EPS0 );
       strainIncrement[1] = 0;
@@ -215,7 +220,7 @@ void TriaxialDriver::runMixedControlTest( SOLID_TYPE & solid, arrayView2d< real6
       integer k = 0;
       integer cuts = 0;
 
-      for(; k<m_maxIter; ++k )
+      for(; k < maxIter; ++k )
       {
         updates.smallStrainUpdate( ei, 0, strainIncrement, stress, stiffness );
 
@@ -226,11 +231,11 @@ void TriaxialDriver::runMixedControlTest( SOLID_TYPE & solid, arrayView2d< real6
           normZero = norm;
         }
 
-        if( norm < m_newtonTol ) // success
+        if( norm < newtonTol ) // success
         {
           break;
         }
-        else if( k > 0 && norm > normZero && cuts < m_maxCuts ) // backtrack by half delta
+        else if( k > 0 && norm > normZero && cuts < maxCuts ) // backtrack by half delta
         {
           cuts++;
           deltaStrainIncrement *= 0.5;
@@ -254,7 +259,7 @@ void TriaxialDriver::runMixedControlTest( SOLID_TYPE & solid, arrayView2d< real6
       table( n, ITER ) = k;
       table( n, NORM ) = norm;
 
-      if( norm > m_newtonTol )
+      if( norm > newtonTol )
       {
         break;
       }
@@ -264,9 +269,13 @@ void TriaxialDriver::runMixedControlTest( SOLID_TYPE & solid, arrayView2d< real6
 
 
 template< typename SOLID_TYPE >
-void TriaxialDriver::runStressControlTest( SOLID_TYPE & solid, arrayView2d< real64 > & table )
+void TriaxialDriver::runStressControlTest( SOLID_TYPE & solid, arrayView2d< real64 > const & table )
 {
   typename SOLID_TYPE::KernelWrapper updates = solid.createKernelUpdates();
+  integer const numSteps = m_numSteps;
+  integer const maxIter = m_maxIter;
+  integer const maxCuts = m_maxCuts;
+  real64 const newtonTol = m_newtonTol;
 
   forAll< parallelDevicePolicy<> >( 1, [=]  GEOSX_HOST_DEVICE ( integer const ei )
   {
@@ -279,13 +288,13 @@ void TriaxialDriver::runStressControlTest( SOLID_TYPE & solid, arrayView2d< real
     real64 jacobian[2][2] = {{}};
 
     real64 scale = 0;
-    for( integer n=1; n<=m_numSteps; ++n )
+    for( integer n = 1; n <= numSteps; ++n )
     {
       scale += fabs( table( n, SIG0 )) + fabs( table( n, SIG1 )) + fabs( table( n, SIG2 ));
     }
-    scale = 3*m_numSteps / scale;
+    scale = 3 * numSteps / scale;
 
-    for( integer n=1; n<=m_numSteps; ++n )
+    for( integer n = 1; n <= numSteps; ++n )
     {
       //   std::cout<<"time step="<<n<<std::endl;
       strainIncrement[0] = 0;
@@ -296,7 +305,7 @@ void TriaxialDriver::runStressControlTest( SOLID_TYPE & solid, arrayView2d< real
       integer k = 0;
       integer cuts = 0;
 
-      for(; k<m_maxIter; ++k )
+      for(; k < maxIter; ++k )
       {
         updates.smallStrainUpdate( ei, 0, strainIncrement, stress, stiffness );
 
@@ -312,11 +321,11 @@ void TriaxialDriver::runStressControlTest( SOLID_TYPE & solid, arrayView2d< real
           normZero = norm;
         }
 
-        if( norm < m_newtonTol ) // success
+        if( norm < newtonTol ) // success
         {
           break;
         }
-        else if( k > 0 && norm > normZero && cuts < m_maxCuts ) // backtrack by half delta
+        else if( k > 0 && norm > normZero && cuts < maxCuts ) // backtrack by half delta
         {
           cuts++;
           deltaStrainIncrement[0] *= 0.5;
@@ -354,7 +363,7 @@ void TriaxialDriver::runStressControlTest( SOLID_TYPE & solid, arrayView2d< real
       table( n, ITER ) = k;
       table( n, NORM ) = norm;
 
-      if( norm > m_newtonTol )
+      if( norm > newtonTol )
       {
         break;
       }

src/coreComponents/constitutive/solid/TriaxialDriver.hpp

@@ -64,23 +64,23 @@ class TriaxialDriver : public TaskBase
    * @param table Table with stress / strain time history
    */
   template< typename SOLID_TYPE >
-  void runStrainControlTest( SOLID_TYPE & solid, arrayView2d< real64 > & table );
+  void runStrainControlTest( SOLID_TYPE & solid, arrayView2d< real64 > const & table );
 
   /**
    * @brief Run a stress-controlled test using loading protocol in table
    * @param solid Solid constitutive model
    * @param table Table with stress / strain time history
    */
   template< typename SOLID_TYPE >
-  void runStressControlTest( SOLID_TYPE & solid, arrayView2d< real64 > & table );
+  void runStressControlTest( SOLID_TYPE & solid, arrayView2d< real64 > const & table );
 
   /**
    * @brief Run a mixed stress/strain-controlled test using loading protocol in table
    * @param solid Solid constitutive model
    * @param table Table with stress / strain time history
    */
   template< typename SOLID_TYPE >
-  void runMixedControlTest( SOLID_TYPE & solid, arrayView2d< real64 > & table );
+  void runMixedControlTest( SOLID_TYPE & solid, arrayView2d< real64 > const & table );
 
   /**
    * @brief Validate results by checking residual and removing erroneous data
@@ -114,7 +114,7 @@ class TriaxialDriver : public TaskBase
     constexpr static char const * baselineString() { return "baseline"; }
   };
 
-  integer m_numSteps;              ///< Number of load steps
+  integer m_numSteps;          ///< Number of load steps
   string m_solidMaterialName;  ///< Material identifier
   string m_mode;               ///< Test mode: strainControl, stressControl, mixedControl
   string m_axialFunctionName;  ///< Time-dependent function controlling axial stress or strain (depends on test mode)

src/coreComponents/dataRepository/BufferOps_inline.hpp

@@ -406,7 +406,7 @@ localIndex Unpack( buffer_unit_type const * & buffer,
 {
   ArrayOfArrays< T > varAsArray;
   localIndex sizeOfUnpackedChars = Unpack( buffer, varAsArray );
-  var.assimilate( std::move( varAsArray ), LvArray::sortedArrayManipulation::SORTED_UNIQUE );
+  var.template assimilate< parallelHostPolicy >( std::move( varAsArray ), LvArray::sortedArrayManipulation::SORTED_UNIQUE );
   return sizeOfUnpackedChars;
 }