Add Tpetra-based Trilinos interface. Refactor LAI vector unit tests.

GEOS-DEV · Aug 8, 2020 · a680573 · a680573
1 parent 5adb70a
commit a680573
Show file tree

Hide file tree

Showing 62 changed files with 4,060 additions and 1,126 deletions.
diff --git a/.travis.yml b/.travis.yml
@@ -2,7 +2,7 @@ language: minimal
 
 env:
   global:
-  - GEOSX_TPL_TAG=119-434
+  - GEOSX_TPL_TAG=120-440
 
 # The integrated test repository contains large data (using git lfs) and we do not use them here.
 # To save time (and money) we do not let travis automatically clone all our (lfs) subrepositories and do it by hand.

diff --git a/src/cmake/GeosxOptions.cmake b/src/cmake/GeosxOptions.cmake
@@ -49,20 +49,22 @@ option( ENABLE_PETSC "Enables PETSC" OFF )
 
 ### LAI SETUP ###
 
-set( supported_LAI Trilinos Hypre Petsc )
+set( supported_LAI Trilinos TrilinosTpetra Hypre Petsc )
 set( GEOSX_LA_INTERFACE "Trilinos" CACHE STRING "Linear algebra interface to use in solvers" )
-message( STATUS "GEOSX_LA_INTERFACE = ${GEOSX_LA_INTERFACE}" )
+message( "GEOSX_LA_INTERFACE = ${GEOSX_LA_INTERFACE}" )
 
 if( NOT ( GEOSX_LA_INTERFACE IN_LIST supported_LAI ) )
   message( FATAL_ERROR "GEOSX_LA_INTERFACE must be one of: ${supported_LAI}" )
 endif()
 
 string( TOUPPER "${GEOSX_LA_INTERFACE}" upper_LAI )
-if( NOT ENABLE_${upper_LAI} )
-  message( FATAL_ERROR "${GEOSX_LA_INTERFACE} LA interface is selected, but ENABLE_${upper_LAI} is OFF" )
-endif()
 option( GEOSX_LA_INTERFACE_${upper_LAI} "${upper_LAI} LA interface is selected" ON )
 
+string( REPLACE "TPETRA" "" dep_LAI "${upper_LAI}" )
+if( NOT ENABLE_${dep_LAI} )
+  message( FATAL_ERROR "${GEOSX_LA_INTERFACE} LA interface is selected, but ENABLE_${dep_LAI} is OFF" )
+endif()
+
 ### MPI/OMP/CUDA SETUP ###
 
 option( ENABLE_MPI "" ON )

diff --git a/src/coreComponents/LvArray b/src/coreComponents/LvArray
diff --git a/src/coreComponents/common/GeosxConfig.hpp.in b/src/coreComponents/common/GeosxConfig.hpp.in
@@ -69,6 +69,8 @@
 #cmakedefine GEOSX_LA_INTERFACE @GEOSX_LA_INTERFACE@
 /// Macro defined when Trilinos interface is selected
 #cmakedefine GEOSX_LA_INTERFACE_TRILINOS
+/// Macro defined when Trilinos interface is selected
+#cmakedefine GEOSX_LA_INTERFACE_TRILINOSTPETRA
 /// Macro defined when Hypre interface is selected
 #cmakedefine GEOSX_LA_INTERFACE_HYPRE
 /// Macro defined when PETSc interface is selected

diff --git a/src/coreComponents/finiteVolume/CMakeLists.txt b/src/coreComponents/finiteVolume/CMakeLists.txt
@@ -26,9 +26,9 @@ set( finiteVolume_sources
    )
 
 if ( BUILD_OBJ_LIBS )
-  set( dependencyList dataRepository codingUtilities managers)
+  set( dependencyList dataRepository codingUtilities managers )
 else()
-  set( dependencyList common)
+  set( dependencyList common managers )
 endif()
 
 if ( ENABLE_OPENMP )

diff --git a/src/coreComponents/linearAlgebra/CMakeLists.txt b/src/coreComponents/linearAlgebra/CMakeLists.txt
@@ -63,16 +63,27 @@ if( ENABLE_TRILINOS )
           interfaces/trilinos/EpetraMatrix.hpp
           interfaces/trilinos/EpetraVector.hpp
           interfaces/trilinos/EpetraUtils.hpp
+          interfaces/trilinos/TpetraMatrix.hpp
+          interfaces/trilinos/TpetraVector.hpp
+          interfaces/trilinos/TpetraUtils.hpp
           interfaces/trilinos/TrilinosPreconditioner.hpp
+          interfaces/trilinos/TrilinosTpetraPreconditioner.hpp
           interfaces/trilinos/TrilinosSolver.hpp
-          interfaces/trilinos/TrilinosInterface.hpp )
+          interfaces/trilinos/TrilinosTpetraSolver.hpp
+          interfaces/trilinos/TrilinosInterface.hpp
+          interfaces/trilinos/TrilinosTpetraInterface.hpp )
 
     list( APPEND linearAlgebra_sources
           interfaces/trilinos/EpetraMatrix.cpp
           interfaces/trilinos/EpetraVector.cpp
+          interfaces/trilinos/TpetraMatrix.cpp
+          interfaces/trilinos/TpetraVector.cpp
           interfaces/trilinos/TrilinosPreconditioner.cpp
+          interfaces/trilinos/TrilinosTpetraPreconditioner.cpp
           interfaces/trilinos/TrilinosSolver.cpp
-          interfaces/trilinos/TrilinosInterface.cpp )
+          interfaces/trilinos/TrilinosTpetraSolver.cpp
+          interfaces/trilinos/TrilinosInterface.cpp
+          interfaces/trilinos/TrilinosTpetraInterface.cpp )
 
     list( APPEND dependencyList trilinos )
 

diff --git a/src/coreComponents/linearAlgebra/DofManager.cpp b/src/coreComponents/linearAlgebra/DofManager.cpp
@@ -1983,6 +1983,7 @@ void DofManager::printFieldInfo( std::ostream & os ) const
 
 #ifdef GEOSX_USE_TRILINOS
 MAKE_DOFMANAGER_METHOD_INST( TrilinosInterface )
+MAKE_DOFMANAGER_METHOD_INST( TrilinosTpetraInterface )
 #endif
 
 #ifdef GEOSX_USE_HYPRE

diff --git a/src/coreComponents/linearAlgebra/interfaces/InterfaceTypes.hpp b/src/coreComponents/linearAlgebra/interfaces/InterfaceTypes.hpp
@@ -23,6 +23,7 @@
 
 #ifdef GEOSX_USE_TRILINOS
 #include "linearAlgebra/interfaces/trilinos/TrilinosInterface.hpp"
+#include "linearAlgebra/interfaces/trilinos/TrilinosTpetraInterface.hpp"
 #endif
 
 #ifdef GEOSX_USE_HYPRE
@@ -62,6 +63,7 @@ inline void setupLAI( int & argc, char * * & argv )
 {
 #ifdef GEOSX_USE_TRILINOS
   TrilinosInterface::initialize( argc, argv );
+  TrilinosTpetraInterface::initialize( argc, argv );
 #endif
 #ifdef GEOSX_USE_HYPRE
   HypreInterface::initialize( argc, argv );
@@ -78,6 +80,7 @@ inline void finalizeLAI()
 {
 #ifdef GEOSX_USE_TRILINOS
   TrilinosInterface::finalize();
+  TrilinosTpetraInterface::finalize();
 #endif
 #ifdef GEOSX_USE_HYPRE
   HypreInterface::finalize();

diff --git a/src/coreComponents/linearAlgebra/interfaces/MatrixBase.hpp b/src/coreComponents/linearAlgebra/interfaces/MatrixBase.hpp
@@ -21,7 +21,6 @@
 
 #include "linearAlgebra/common.hpp"
 #include "linearAlgebra/interfaces/LinearOperator.hpp"
-//#include "LvArray/src/streamIO.hpp"
 
 namespace geosx
 {
@@ -438,36 +437,6 @@ class MatrixBase : public virtual LinearOperator< VECTOR >
                        arraySlice1d< globalIndex const > const & colIndices,
                        arraySlice2d< real64 const, MatrixLayout::ROW_MAJOR > const & values ) = 0;
 
-  /**
-   * @brief Add a dense block of values.
-   * @param rowIndices Global row indices
-   * @param colIndices Global col indices
-   * @param values Dense local matrix of values
-   */
-  virtual void add( arraySlice1d< globalIndex const > const & rowIndices,
-                    arraySlice1d< globalIndex const > const & colIndices,
-                    arraySlice2d< real64 const, MatrixLayout::COL_MAJOR > const & values ) = 0;
-
-  /**
-   * @brief Set a dense block of values.
-   * @param rowIndices Global row indices
-   * @param colIndices Global col indices
-   * @param values Dense local matrix of values
-   */
-  virtual void set( arraySlice1d< globalIndex const > const & rowIndices,
-                    arraySlice1d< globalIndex const > const & colIndices,
-                    arraySlice2d< real64 const, MatrixLayout::COL_MAJOR > const & values ) = 0;
-
-  /**
-   * @brief Insert a dense block of values.
-   * @param rowIndices Global row indices
-   * @param colIndices Global col indices
-   * @param values Dense local matrix of values
-   */
-  virtual void insert( arraySlice1d< globalIndex const > const & rowIndices,
-                       arraySlice1d< globalIndex const > const & colIndices,
-                       arraySlice2d< real64 const, MatrixLayout::COL_MAJOR > const & values ) = 0;
-
   /**
    * @brief Add a dense block of values.
    * @param rowIndices Global row indices

diff --git a/src/coreComponents/linearAlgebra/interfaces/VectorBase.hpp b/src/coreComponents/linearAlgebra/interfaces/VectorBase.hpp
@@ -149,7 +149,7 @@ class VectorBase
    * @param localValues local data to put into vector
    * @param comm MPI communicator to use
    */
-  virtual void create( arrayView1d< real64 const > const & localValues, MPI_Comm const & comm ) = 0;
+  virtual void createWithLocalValues( arrayView1d< real64 > const & localValues, MPI_Comm const & comm ) = 0;
 
   ///@}
 
@@ -414,6 +414,7 @@ class VectorBase
    */
   virtual void extract( arrayView1d< real64 > const & localVector ) const
   {
+    GEOSX_LAI_ASSERT_EQ( localSize(), localVector.size() );
     real64 const * const data = extractLocalVector();
     forAll< parallelHostPolicy >( localSize(), [=] ( localIndex const k )
     {

diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.cpp
@@ -25,6 +25,7 @@
 #include "HypreUtils.hpp"
 
 #include <iomanip>
+#include <numeric>
 
 namespace geosx
 {
@@ -444,53 +445,6 @@ void HypreMatrix::insert( arraySlice1d< globalIndex const > const & rowIndices,
   }
 }
 
-namespace
-{
-
-template< typename T, int SRC_USD, int DST_USD >
-static void convertArrayLayout( arraySlice2d< T const, SRC_USD > const & src,
-                                arraySlice2d< T, DST_USD > const & dst )
-{
-  GEOSX_ASSERT( src.size( 0 ) == dst.size( 0 ) );
-  GEOSX_ASSERT( src.size( 1 ) == dst.size( 1 ) );
-  for( localIndex i = 0; i < src.size( 0 ); ++i )
-  {
-    for( localIndex j = 0; j < src.size( 1 ); ++j )
-    {
-      dst( i, j ) = src( i, j );
-    }
-  }
-}
-
-}
-
-void HypreMatrix::add( arraySlice1d< globalIndex const > const & rowIndices,
-                       arraySlice1d< globalIndex const > const & colIndices,
-                       arraySlice2d< real64 const, MatrixLayout::COL_MAJOR > const & values )
-{
-  array2d< real64, MatrixLayout::ROW_MAJOR_PERM > valuesRowMajor( rowIndices.size(), colIndices.size() );
-  convertArrayLayout( values, valuesRowMajor.toSlice() );
-  add( rowIndices, colIndices, valuesRowMajor );
-}
-
-void HypreMatrix::set( arraySlice1d< globalIndex const > const & rowIndices,
-                       arraySlice1d< globalIndex const > const & colIndices,
-                       arraySlice2d< real64 const, MatrixLayout::COL_MAJOR > const & values )
-{
-  array2d< real64, MatrixLayout::ROW_MAJOR_PERM > valuesRowMajor( rowIndices.size(), colIndices.size() );
-  convertArrayLayout( values, valuesRowMajor.toSlice() );
-  set( rowIndices, colIndices, valuesRowMajor );
-}
-
-void HypreMatrix::insert( arraySlice1d< globalIndex const > const & rowIndices,
-                          arraySlice1d< globalIndex const > const & colIndices,
-                          arraySlice2d< real64 const, MatrixLayout::COL_MAJOR > const & values )
-{
-  array2d< real64, MatrixLayout::ROW_MAJOR_PERM > valuesRowMajor( rowIndices.size(), colIndices.size() );
-  convertArrayLayout( values, valuesRowMajor.toSlice() );
-  insert( rowIndices, colIndices, valuesRowMajor );
-}
-
 void HypreMatrix::add( globalIndex const * rowIndices,
                        globalIndex const * colIndices,
                        real64 const * values,
@@ -905,10 +859,7 @@ localIndex HypreMatrix::maxRowLength() const
   array1d< HYPRE_BigInt > rows( nrows );
   array1d< HYPRE_Int > ncols( nrows );
 
-  for( HYPRE_Int i = 0; i < nrows; ++i )
-  {
-    rows[i] = ilower() + LvArray::integerConversion< HYPRE_BigInt >( i );
-  }
+  std::iota( rows.begin(), rows.end(), ilower() );
 
   GEOSX_LAI_CHECK_ERROR( HYPRE_IJMatrixGetRowCounts( m_ij_mat,
                                                      nrows,

diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.hpp
@@ -177,18 +177,6 @@ class HypreMatrix final : public virtual LinearOperator< HypreVector >,
                        arraySlice1d< globalIndex const > const & colIndices,
                        arraySlice2d< real64 const, MatrixLayout::ROW_MAJOR > const & values ) override;
 
-  virtual void add( arraySlice1d< globalIndex const > const & rowIndices,
-                    arraySlice1d< globalIndex const > const & colIndices,
-                    arraySlice2d< real64 const, MatrixLayout::COL_MAJOR > const & values ) override;
-
-  virtual void set( arraySlice1d< globalIndex const > const & rowIndices,
-                    arraySlice1d< globalIndex const > const & colIndices,
-                    arraySlice2d< real64 const, MatrixLayout::COL_MAJOR > const & values ) override;
-
-  virtual void insert( arraySlice1d< globalIndex const > const & rowIndices,
-                       arraySlice1d< globalIndex const > const & colIndices,
-                       arraySlice2d< real64 const, MatrixLayout::COL_MAJOR > const & values ) override;
-
   virtual void add( globalIndex const * rowIndices,
                     globalIndex const * colIndices,
                     real64 const * values,

diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.cpp
@@ -22,7 +22,6 @@
 #include "linearAlgebra/DofManager.hpp"
 #include "linearAlgebra/interfaces/hypre/HypreUtils.hpp"
 #include "linearAlgebra/utilities/LinearSolverParameters.hpp"
-#include "linearAlgebra/utilities/LAIHelperFunctions.hpp"
 
 #include <_hypre_utilities.h>
 #include <_hypre_parcsr_ls.h>
@@ -67,6 +66,10 @@ HyprePreconditioner::HyprePreconditioner( LinearSolverParameters params,
     {
       createILUT();
     }
+    else if( m_parameters.preconditionerType == "direct" )
+    {
+      createDirect();
+    }
     else
     {
       GEOSX_ERROR( "Unsupported preconditioner type: " << m_parameters.preconditionerType );
@@ -534,6 +537,42 @@ void HyprePreconditioner::createILUT()
   m_functions->destroy = HYPRE_ILUDestroy;
 }
 
+namespace internal
+{
+
+static HYPRE_Int
+HYPRE_SLUDistSetup( HYPRE_Solver solver,
+                    HYPRE_ParCSRMatrix A,
+                    HYPRE_ParVector GEOSX_UNUSED_PARAM( b ),
+                    HYPRE_ParVector GEOSX_UNUSED_PARAM( x ) )
+{
+  return hypre_SLUDistSetup( &solver, A, 0 );
+}
+
+static HYPRE_Int
+HYPRE_SLUDistSolve( HYPRE_Solver solver,
+                    HYPRE_ParCSRMatrix GEOSX_UNUSED_PARAM( A ),
+                    HYPRE_ParVector b,
+                    HYPRE_ParVector x )
+{
+  return hypre_SLUDistSolve( solver, b, x );
+}
+
+static HYPRE_Int
+HYPRE_SLUDistDestroy( HYPRE_Solver solver )
+{
+  return hypre_SLUDistDestroy( solver );
+}
+
+}
+
+void HyprePreconditioner::createDirect()
+{
+  m_functions->setup = internal::HYPRE_SLUDistSetup;
+  m_functions->apply = internal::HYPRE_SLUDistSolve;
+  m_functions->destroy = internal::HYPRE_SLUDistDestroy;
+}
+
 void HyprePreconditioner::compute( Matrix const & mat )
 {
   PreconditionerBase::compute( mat );

diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.hpp
@@ -123,6 +123,8 @@ class HyprePreconditioner final : public PreconditionerBase< HypreInterface >
 
   void createILUT();
 
+  void createDirect();
+
   /// Parameters for all preconditioners
   LinearSolverParameters m_parameters;
 
@@ -135,7 +137,7 @@ class HyprePreconditioner final : public PreconditionerBase< HypreInterface >
   /// Pointers to hypre functions to setup/solve/destroy preconditioner
   std::unique_ptr< HyprePrecFuncs > m_functions;
 
-  // Pointer to preconditioner auxiliary data
+  /// Pointer to preconditioner auxiliary data
   std::unique_ptr< HyprePrecAuxData > m_auxData;
 };
 

diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreVector.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreVector.cpp
@@ -157,8 +157,8 @@ void HypreVector::createWithGlobalSize( globalIndex const globalSize,
   finalize( m_ij_vector, m_par_vector );
 }
 
-void HypreVector::create( arrayView1d< real64 const > const & localValues,
-                          MPI_Comm const & comm )
+void HypreVector::createWithLocalValues( arrayView1d< real64 > const & localValues,
+                                         MPI_Comm const & comm )
 {
   GEOSX_LAI_ASSERT( closed() );
 

diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreVector.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreVector.hpp
@@ -117,8 +117,8 @@ class HypreVector final : private VectorBase< HypreVector >
   virtual void createWithGlobalSize( globalIndex const globalSize,
                                      MPI_Comm const & comm ) override;
 
-  virtual void create( arrayView1d< real64 const > const & localValues,
-                       MPI_Comm const & comm ) override;
+  virtual void createWithLocalValues( arrayView1d< real64 > const & localValues,
+                                      MPI_Comm const & comm ) override;
 
   virtual void open() override;