Remove R1Tensor (#1180)

* Removed Cross products. Fixed BoundedPlane, ThickPlane, ComputationalGeometry. Started removal of Dot from solvers and surfGenerators. * Clean up more R1Tensor uses. Get rid or realT. * A few more cleanups. * Remove most R1Tensor uses in flow solvers * surfaceGen r1Tensor operators * R1Tensor initialization in surfGen * removed R1Tensor use from all remaining solvers. * WIP change R1Tensor in EmbSurfSubregion * Templated functions on array_type in CompGeometry. * Add Simple Tensor type for input * Bugfix local var initialization Co-authored-by: Sergey Klevtsov <[email protected]> Co-authored-by: Randolph R. Settgast <[email protected]>
GEOS-DEV · Nov 26, 2020 · 1e1a3e6 · 1e1a3e6
1 parent ecaa0f2
commit 1e1a3e6
Show file tree

Hide file tree

Showing 103 changed files with 1,853 additions and 3,679 deletions.
diff --git a/.gitmodules b/.gitmodules
@@ -1,22 +1,22 @@
 [submodule "src/coreComponents/LvArray"]
 	path = src/coreComponents/LvArray
-	url = https://github.com/GEOSX/LvArray.git
+	url = ../LvArray.git
 [submodule "src/cmake/blt"]
 	path = src/cmake/blt
-	url = https://github.com/LLNL/blt.git
+	url = ../../LLNL/blt.git
 [submodule "src/externalComponents/PVTPackage"]
 	path = src/externalComponents/PVTPackage
-	url = https://github.com/GEOSX/PVTPackage.git
+	url = ../PVTPackage.git
 [submodule "integratedTests"]
 	path = integratedTests
-	url = [email protected]:GEOSX/integratedTests.git
+	url = ../integratedTests.git
 [submodule "src/externalComponents/PAMELA"]
 	path = src/externalComponents/PAMELA
-	url = https://github.com/GEOSX/PAMELA.git
+	url = ../PAMELA.git
 	shallow = true
 [submodule "src/coreComponents/fileIO/coupling/hdf5_interface"]
 	path = src/coreComponents/fileIO/coupling/hdf5_interface
-	url = https://github.com/GEOSX/hdf5_interface.git
+	url = ../hdf5_interface.git
 [submodule "src/coreComponents/physicsSolvers/GEOSX_PTP"]
 	path = src/coreComponents/physicsSolvers/GEOSX_PTP
-	url = https://github.com/GEOSX/GEOSX_PTP.git
+	url = ../GEOSX_PTP.git
diff --git a/src/coreComponents/LvArray b/src/coreComponents/LvArray
diff --git a/src/coreComponents/codingUtilities/SFINAE_Macros.hpp b/src/coreComponents/codingUtilities/SFINAE_Macros.hpp
@@ -13,11 +13,11 @@
  */
 
 /**
- * @file sfinae.hpp
+ * @file SFINAE_Macros.hpp
  */
 
-#ifndef SRC_CODINGUTILITIES_SFINAE_HPP_
-#define SRC_CODINGUTILITIES_SFINAE_HPP_
+#ifndef GEOSX_CODINGUTILITIES_SFINAE_MACROS_HPP_
+#define GEOSX_CODINGUTILITIES_SFINAE_MACROS_HPP_
 
 #include "LvArray/src/Macros.hpp"
 #include "LvArray/src/typeManipulation.hpp"
@@ -50,11 +50,13 @@
   IS_VALID_EXPRESSION( HasMemberFunction_ ## NAME, CLASS, \
                        std::is_convertible< decltype( std::declval< CLASS >().NAME( __VA_ARGS__ ) ), RTYPE >::value )
 
+
 /**
  * @brief Macro that expands to a static constexpr bool templated on a type that is only true when
  *        the type has a an alias @p NAME. The name of the boolean variable is HasAlias_ ## @p NAME.
  */
 #define HAS_ALIAS( NAME ) \
   IS_VALID_EXPRESSION( HasAlias_ ## NAME, CLASS, !std::is_enum< typename CLASS::NAME >::value )
 
-#endif /* SRC_CODINGUTILITIES_SFINAE_HPP_ */
+
+#endif /* GEOSX_CODINGUTILITIES_SFINAE_MACROS_HPP_ */
diff --git a/src/coreComponents/codingUtilities/Utilities.hpp b/src/coreComponents/codingUtilities/Utilities.hpp
@@ -20,22 +20,11 @@
 #define GEOSX_CODINGUTILITIES_UTILITIES_H_
 
 #include "common/DataTypes.hpp"
+#include "LvArray/src/limits.hpp"
 
 namespace geosx
 {
 
-/**
- * @brief GPU-friendly analogue of std::numeric_limits
- * @tparam T type of numeric value
- */
-template< typename T >
-struct NumericTraits
-{
-  static constexpr T min = std::numeric_limits< T >::lowest();
-  static constexpr T max = std::numeric_limits< T >::max();
-  static constexpr T eps = std::numeric_limits< T >::epsilon();
-};
-
 /**
  * @brief Compare two real values with a tolerance.
  * @tparam type of real value
@@ -48,7 +37,7 @@ template< typename T >
 GEOSX_FORCE_INLINE GEOSX_HOST_DEVICE constexpr
 bool isEqual( T const val1, T const val2, T const relTol = 0.0 )
 {
-  T const absTol = ( relTol > 10 * NumericTraits< T >::eps ) ? relTol * ( fabs( val1 ) + fabs( val2 ) ) * 0.5 : 0.0;
+  T const absTol = ( relTol > 10 * LvArray::NumericLimits< T >::epsilon ) ? relTol * ( fabs( val1 ) + fabs( val2 ) ) * 0.5 : 0.0;
   return ( val2 - absTol ) <= val1 && val1 <= ( val2 + absTol );
 }
 
@@ -61,7 +50,7 @@ bool isEqual( T const val1, T const val2, T const relTol = 0.0 )
  */
 template< typename T >
 GEOSX_FORCE_INLINE GEOSX_HOST_DEVICE constexpr
-bool isZero( T const val, T const tol = NumericTraits< T >::eps )
+bool isZero( T const val, T const tol = LvArray::NumericLimits< T >::epsilon )
 {
   return -tol <= val && val <= tol;
 }

diff --git a/src/coreComponents/codingUtilities/tests/testGeosxTraits.cpp b/src/coreComponents/codingUtilities/tests/testGeosxTraits.cpp
@@ -28,7 +28,6 @@ TEST( testGeosxTraits, HasOperatorPlusEquals )
 {
   static_assert( HasOperatorPlusEquals< int >, "Should be true." );
   static_assert( HasOperatorPlusEquals< std::string >, "Should be true." );
-  static_assert( HasOperatorPlusEquals< R1Tensor >, "Should be true." );
 
   static_assert( !HasOperatorPlusEquals< int const >, "Should be false." );
   static_assert( !HasOperatorPlusEquals< std::vector< int > >, "Should be false." );
@@ -39,8 +38,6 @@ TEST( testGeosxTraits, HasOperatorPlusEquals2 )
 {
   static_assert( HasOperatorPlusEquals2< double, int >, "Should be true." );
   static_assert( HasOperatorPlusEquals2< double, int const >, "Should be true." );
-  static_assert( HasOperatorPlusEquals2< R1Tensor, int >, "Should be true." );
-  static_assert( HasOperatorPlusEquals2< R1Tensor, R1Tensor >, "Should be true." );
   static_assert( HasOperatorPlusEquals2< std::string, std::string >, "Should be true." );
 
   static_assert( !HasOperatorPlusEquals2< double const, int >, "Should be false." );

diff --git a/src/coreComponents/codingUtilities/traits.hpp b/src/coreComponents/codingUtilities/traits.hpp
@@ -101,7 +101,7 @@ static constexpr bool HasMemberFunction_toView = LvArray::typeManipulation::HasM
 IS_VALID_EXPRESSION_2( CanStreamInto, SRC, DST, std::declval< SRC & >() >> std::declval< DST & >() );
 
 /**
- * @brief Defines a static constexpr bool HasAlias_value_type< @p CLASS >
+ * @brief Defines a static constexpr bool HasMemberType_value_type< @p CLASS >
  *        that is true iff @p CLASS ::value_type is valid and not an enum.
  * @tparam CLASS The type to test.
  */

diff --git a/src/coreComponents/common/CMakeLists.txt b/src/coreComponents/common/CMakeLists.txt
@@ -13,6 +13,7 @@ set(common_headers
     TimingMacros.hpp
     Logger.hpp
     DataLayouts.hpp
+    Tensor.hpp
    )
 
 #

diff --git a/src/coreComponents/common/DataTypes.hpp b/src/coreComponents/common/DataTypes.hpp
@@ -27,7 +27,8 @@
 #include "common/GeosxMacros.hpp"
 #include "common/BufferAllocator.hpp"
 #include "common/DataLayouts.hpp"
-#include "Logger.hpp"
+#include "common/Tensor.hpp"
+#include "common/Logger.hpp"
 #include "LvArray/src/Macros.hpp"
 #include "LvArray/src/Array.hpp"
 #include "LvArray/src/ArrayOfArrays.hpp"
@@ -39,7 +40,6 @@
 #include "LvArray/src/StackBuffer.hpp"
 #include "LvArray/src/ChaiBuffer.hpp"
 
-#include "math/TensorT/TensorT.h"
 #include "Path.hpp"
 
 // TPL includes
@@ -189,6 +189,9 @@ using StackArray = LvArray::StackArray< T, NDIM, PERMUTATION, localIndex, MAXSIZ
  */
 ///@{
 
+/// Alias for a local (stack-based) rank-1 tensor type
+using R1Tensor = Tensor< real64, 3 >;
+
 /// Alias for 1D array.
 template< typename T >
 using array1d = Array< T, 1 >;

diff --git a/src/coreComponents/common/Logger.hpp b/src/coreComponents/common/Logger.hpp
@@ -66,6 +66,9 @@
  * @param EXP an expression that will be evaluated as a predicate
  * @param msg a message to log (any expression that can be stream inserted)
  */
+#if defined(__CUDA_ARCH__)
+#define GEOSX_LOG_RANK_IF( EXP, msg )
+#else
 #define GEOSX_LOG_RANK_IF( EXP, msg ) \
   do { \
     if( EXP ) \
@@ -75,6 +78,7 @@
       *logger::internal::rankStream << oss.str() << std::endl; \
     } \
   } while( false )
+#endif
 
 /**
  * @brief Log a message to the rank output stream.
@@ -93,7 +97,11 @@
  * @param EXP an expression that will be evaluated as a predicate
  * @param msg a message to log (any expression that can be stream inserted)
  */
+#if defined(__CUDA_ARCH__)
+#define GEOSX_ERROR_IF( EXP, msg ) LVARRAY_ERROR_IF( EXP, msg )
+#else
 #define GEOSX_ERROR_IF( EXP, msg ) LVARRAY_ERROR_IF( EXP, "***** Rank " << ::geosx::logger::internal::rankString << ": " << msg )
+#endif
 
 /**
  * @brief Raise a hard error and terminate the program.

diff --git a/src/coreComponents/common/Tensor.hpp b/src/coreComponents/common/Tensor.hpp
@@ -0,0 +1,148 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2018-2020 Total, S.A
+ * Copyright (c) 2019-     GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+#ifndef GEOSX_COMMON_TENSOR_HPP_
+#define GEOSX_COMMON_TENSOR_HPP_
+
+#include "common/GeosxMacros.hpp"
+
+namespace geosx
+{
+
+/**
+ * @brief Lightweight wrapper around a c-array
+ * @tparam T data type
+ * @tparam SIZE_TPARAM number of values
+ */
+template< typename T, int SIZE_TPARAM >
+class Tensor
+{
+public:
+
+  static_assert( SIZE_TPARAM > 0, "Tensor size must be a positive value" );
+
+  /// Alias for type template parameter
+  using value_type = T;
+
+  /// Alias for size template parameter
+  static constexpr int SIZE = SIZE_TPARAM;
+
+  /**
+   * @brief Const element access.
+   * @param i element index
+   * @return const reference to the i-th element
+   */
+  GEOSX_HOST_DEVICE
+  GEOSX_FORCE_INLINE
+  T const & operator[]( std::ptrdiff_t const i ) const
+  {
+    return data[i];
+  }
+
+  /**
+   * @brief Non-const element access.
+   * @param i element index
+   * @return reference to the i-th element
+   */
+  GEOSX_HOST_DEVICE
+  GEOSX_FORCE_INLINE
+  T & operator[]( std::ptrdiff_t const i )
+  {
+    return data[i];
+  }
+
+  /**
+   * @brief Equality comparison operator.
+   * @tparam U dummy parameter to enable SFINAE, do not provide explicitly
+   * @param rhs the right-hand side value to compare to
+   * @return true iff @code data[i] == rhs.data[i] @endcode
+   * @note Version for floating point types that avoids direct equality comparison.
+   */
+  template< typename U = T >
+  GEOSX_HOST_DEVICE
+  GEOSX_FORCE_INLINE
+  std::enable_if_t< std::is_floating_point< U >::value, bool >
+  operator==( Tensor< U, SIZE > const & rhs ) const
+  {
+    for( int i = 0; i < SIZE; ++i )
+    {
+      if( (data[i] > rhs.data[i]) || (data[i] < rhs.data[i]) )
+      {
+        return false;
+      }
+    }
+    return true;
+  }
+
+  /**
+   * @brief Equality comparison operator.
+   * @tparam U dummy parameter to enable SFINAE, do not provide explicitly
+   * @param rhs the right-hand side value to compare to
+   * @return true iff @code data[i] == rhs.data[i] @endcode
+   * @note Version for all types except floating point.
+   */
+  template< typename U = T >
+  GEOSX_HOST_DEVICE
+  GEOSX_FORCE_INLINE
+  std::enable_if_t< !std::is_floating_point< U >::value, bool >
+  operator==( Tensor< U, SIZE > const & rhs ) const
+  {
+    for( int i = 0; i < SIZE; ++i )
+    {
+      if( data[i] != rhs.data[i] )
+      {
+        return false;
+      }
+    }
+    return true;
+  }
+
+  /**
+   * @brief Returns the size of the tensor
+   * @param junk Unused
+   * @return The value of the template parameter SIZE.
+   */
+  GEOSX_HOST_DEVICE
+  GEOSX_FORCE_INLINE
+  constexpr int size( int junk ) const
+  {
+    GEOSX_UNUSED_VAR( junk )
+    return SIZE;
+  }
+
+  /// Underlying array
+  T data[SIZE] = {};
+
+private:
+
+  /**
+   * @brief Stream insertion operator for Tensor.
+   * @param os the output stream
+   * @param t the tensor value
+   * @return reference to @p os
+   */
+  friend inline std::ostream & operator<<( std::ostream & os, Tensor< T, SIZE > const & t )
+  {
+    os << t.data[0];
+    for( int i = 1; i < SIZE; ++i )
+    {
+      os << ',' << t.data[i];
+    }
+    return os;
+  }
+};
+
+}
+
+#endif //GEOSX_COMMON_TENSOR_HPP_
diff --git a/src/coreComponents/constitutive/fluid/MultiPhaseMultiComponentFluid.cpp b/src/coreComponents/constitutive/fluid/MultiPhaseMultiComponentFluid.cpp
@@ -273,7 +273,7 @@ void MultiPhaseMultiComponentFluidUpdate::Compute( real64 pressure,
       X[ic].m_var = composition[ic];
       X[ic].m_der[ic+1] = 1.0;
 
-      realT const mwInv = 1.0 / m_componentMolarWeight[ic];
+      real64 const mwInv = 1.0 / m_componentMolarWeight[ic];
       C[ic] = X[ic] * mwInv; // this is molality (units of mole/mass)
       totalMolality += C[ic];
     }
@@ -342,7 +342,7 @@ void MultiPhaseMultiComponentFluidUpdate::Compute( real64 pressure,
       for( localIndex ic = 0; ic < NC; ++ic )
       {
 
-        realT compMW = m_componentMolarWeight[ic];
+        real64 compMW = m_componentMolarWeight[ic];
 
         phaseCompFractionTemp[ip][ic] = phaseCompFractionTemp[ip][ic] * compMW /  phaseMW[ip];
 

diff --git a/...coreComponents/constitutive/relativePermeability/BrooksCoreyBakerRelativePermeability.hpp b/...coreComponents/constitutive/relativePermeability/BrooksCoreyBakerRelativePermeability.hpp
@@ -362,13 +362,13 @@ BrooksCoreyBakerRelativePermeabilityUpdate::
   integer const ip_gas   = phaseOrder[RelativePermeabilityBase::PhaseType::GAS];
 
   // if water phase is immobile, then use the two-phase gas-oil data only
-  if( shiftedWaterVolFrac < NumericTraits< real64 >::eps )
+  if( shiftedWaterVolFrac < LvArray::NumericLimits< real64 >::epsilon )
   {
     threePhaseRelPerm = relPerm_go;
     dThreePhaseRelPerm_dVolFrac[ip_oil] = dRelPerm_go_dOilVolFrac;
   }
   // if gas phase is immobile, then use the two-phase water-oil data only
-  else if( gasVolFrac < NumericTraits< real64 >::eps )
+  else if( gasVolFrac < LvArray::NumericLimits< real64 >::epsilon )
   {
     threePhaseRelPerm = relPerm_wo;
     dThreePhaseRelPerm_dVolFrac[ip_oil] = dRelPerm_wo_dOilVolFrac;

diff --git a/...oreComponents/constitutive/relativePermeability/VanGenuchtenBakerRelativePermeability.hpp b/...oreComponents/constitutive/relativePermeability/VanGenuchtenBakerRelativePermeability.hpp
@@ -373,13 +373,13 @@ VanGenuchtenBakerRelativePermeabilityUpdate::
   integer const ip_gas   = phaseOrder[RelativePermeabilityBase::PhaseType::GAS];
 
   // if water phase is immobile, then use the two-phase gas-oil data only
-  if( shiftedWaterVolFrac < NumericTraits< real64 >::eps )
+  if( shiftedWaterVolFrac < LvArray::NumericLimits< real64 >::epsilon )
   {
     threePhaseRelPerm = relPerm_go;
     dThreePhaseRelPerm_dVolFrac[ip_oil] = dRelPerm_go_dOilVolFrac;
   }
   // if gas phase is immobile, then use the two-phase water-oil data only
-  else if( gasVolFrac < NumericTraits< real64 >::eps )
+  else if( gasVolFrac < LvArray::NumericLimits< real64 >::epsilon )
   {
     threePhaseRelPerm = relPerm_wo;
     dThreePhaseRelPerm_dVolFrac[ip_oil] = dRelPerm_wo_dOilVolFrac;
+77 −1		docs/sphinx/developmentAids.rst
+28 −17		scripts/gdb-printers.py
+68 −17		src/genericTensorOps.hpp
+18 −0		src/typeManipulation.hpp
+83 −0		unitTests/testTensorOpsOneSize.cpp