diff --git a/.travis.yml b/.travis.yml index f9a0a530db6..0b5fadc9e18 100644 --- a/.travis.yml +++ b/.travis.yml @@ -2,7 +2,7 @@ language: minimal env: global: - - GEOSX_TPL_TAG=123-507 + - GEOSX_TPL_TAG=120-511 - secure: CGs2uH6efq1Me6xJWRr0BnwtwxoujzlowC4FHXHdWbNOkPsXf7nCgdaW5vthfD3bhnOeEUQSrfxdhTRtyU/NfcKLmKgGBnZOdUG4/JJK4gDSJ2Wp8LZ/mB0QEoODKVxbh+YtoAiHe3y4M9PGCs+wkNDw/3eEU00cK12DZ6gad0RbLjI3xkhEr/ZEZDZkcYg9yHAhl5bmpqoh/6QGnIg8mxIqdAtGDw+6tT0EgUqjeqc5bG5WwsamKzJItHSXD5zx8IJAlgDk4EzEGjZe0m56YnNfb9iwqqUsmL3Cuwgs7ByVDYw78JC5Kv42YqoxA5BxMT2mFsEe37TpYNXlzofU7ma2Duw9DGXWQd4IkTCcBxlyR0I0bfo0TmgO+y7PYG9lIyHPUkENemdozsZcWamqqkqegiEdRhDVYlSRo3mu7iCwTS6ZTALliVyEYjYxYb7oAnR3cNywXjblTCI8oKfgLSY+8WijM9SRl57JruIHLkLMCjmRI+cZBfv5tS2tYQTBPkygGrigrrN77ZiC7/TGyfggSN0+y0oYtOAgqEnBcKcreiibMW7tKcV2Z1RFD9ZvIkSc1EXLUPDP8FX1oyhmqBMqVo8LksrYLDJHQ05+F3YNgl2taSt7uMjQ4e8iZ3/IjFeMnbylDw+cj/RbS520HXsFPbWFm2Pb9pceA9n6GnY= # The integrated test repository contains large data (using git lfs) and we do not use them here. @@ -17,7 +17,7 @@ geosx_before_script: &geosx_before_script - git submodule update --init src/externalComponents/PAMELA - git submodule update --init --recursive src/coreComponents/fileIO/coupling/hdf5_interface - git submodule update --init --recursive src/coreComponents/physicsSolvers/GEOSX_PTP - + geosx_linux_build: &geosx_linux_build services: docker diff --git a/integratedTests b/integratedTests index b28a9fddef3..a6d37a2a07d 160000 --- a/integratedTests +++ b/integratedTests @@ -1 +1 @@ -Subproject commit b28a9fddef3ab3d5238f2599c5d5f895d669e9b1 +Subproject commit a6d37a2a07d7429aad699ed03183d205e30b81c1 diff --git a/src/cmake/GeosxOptions.cmake b/src/cmake/GeosxOptions.cmake index 8bf3e809407..b38a9102cbb 100644 --- a/src/cmake/GeosxOptions.cmake +++ b/src/cmake/GeosxOptions.cmake @@ -51,20 +51,22 @@ option( ENABLE_SUITESPARSE "Enables SUITESPARSE" ON ) ### 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 index f634198bf60..95db802921c 160000 --- a/src/coreComponents/LvArray +++ b/src/coreComponents/LvArray @@ -1 +1 @@ -Subproject commit f634198bf6080ff45198e2b3d263027e2fcd671d +Subproject commit 95db802921cc2da0e43b92763ffecd6a2615e749 diff --git a/src/coreComponents/common/GeosxConfig.hpp.in b/src/coreComponents/common/GeosxConfig.hpp.in index 27ae279a3c5..606ed83ddfc 100644 --- a/src/coreComponents/common/GeosxConfig.hpp.in +++ b/src/coreComponents/common/GeosxConfig.hpp.in @@ -75,6 +75,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/fileIO/schema/docs/LinearSolverParameters.rst b/src/coreComponents/fileIO/schema/docs/LinearSolverParameters.rst index 79c116c2b01..183d28598f1 100644 --- a/src/coreComponents/fileIO/schema/docs/LinearSolverParameters.rst +++ b/src/coreComponents/fileIO/schema/docs/LinearSolverParameters.rst @@ -1,60 +1,86 @@ -=================== =============================================== =========== ======================================================================================================================================================================================================================================================================================================================= -Name Type Default Description -=================== =============================================== =========== ======================================================================================================================================================================================================================================================================================================================= -amgCoarseSolver string direct | AMG coarsest level solver/smoother type - | Available options are: jacobi, gaussSeidel, blockGaussSeidel, chebyshev, direct -amgNumSweeps integer 2 AMG smoother sweeps -amgSmootherType string gaussSeidel | AMG smoother type - | Available options are: jacobi, blockJacobi, gaussSeidel, blockGaussSeidel, chebyshev, icc, ilu, ilut -amgThreshold real64 0 AMG strength-of-connection threshold -directCheckResTol real64 1e-12 Tolerance used to check a direct solver solution -directColPerm geosx_LinearSolverParameters_Direct_ColPerm metis | How to permute the columns. Available options are: - | * none - | * MMD_AtplusA - | * MMD_AtA - | * colAMD - | * metis - | * parmetis -directEquil integer 1 Whether to scale the rows and columns of the matrix -directIterRef integer 1 Whether to perform iterative refinement -directParallel integer 1 Whether to use a parallel solver (instead of a serial one) -directReplTinyPivot integer 1 Whether to replace tiny pivots by sqrt(epsilon)*norm(A) -directRowPerm geosx_LinearSolverParameters_Direct_RowPerm mc64 | How to permute the rows. Available options are: - | * none - | * mc64 -iluFill integer 0 ILU(K) fill factor -iluThreshold real64 0 ILU(T) threshold factor -krylovAdaptiveTol integer 0 Use Eisenstat-Walker adaptive linear tolerance -krylovMaxIter integer 200 Maximum iterations allowed for an iterative solver -krylovMaxRestart integer 200 Maximum iterations before restart (GMRES only) -krylovTol real64 1e-06 | Relative convergence tolerance of the iterative method - | If the method converges, the iterative solution :math:`\mathsf{x}_k` is such that - | the relative residual norm satisfies: - | :math:`\left\lVert \mathsf{b} - \mathsf{A} \mathsf{x}_k \right\rVert_2` < ``krylovTol`` * :math:`\left\lVert\mathsf{b}\right\rVert_2` -krylovWeakestTol real64 0.001 Weakest-allowed tolerance for adaptive method -logLevel integer 0 Log level -preconditionerType geosx_LinearSolverParameters_PreconditionerType iluk | Preconditioner type. Available options are: - | * none - | * jacobi - | * gs - | * sgs - | * iluk - | * ilut - | * icc - | * ict - | * amg - | * mgr - | * block -solverType geosx_LinearSolverParameters_SolverType direct | Linear solver type. Available options are: - | * direct - | * cg - | * gmres - | * fgmres - | * bicgstab - | * preconditioner -stopIfError integer 1 Whether to stop the simulation if the linear solver reports an error -=================== =============================================== =========== ======================================================================================================================================================================================================================================================================================================================= +=================== =============================================== ======= ======================================================================================================================================================================================================================================================================================================================= +Name Type Default Description +=================== =============================================== ======= ======================================================================================================================================================================================================================================================================================================================= +amgCoarseSolver geosx_LinearSolverParameters_PreconditionerType direct | AMG coarsest level solver/smoother type. Valid options (not all may be supported by linear algebra package): + | * none + | * jacobi + | * gs + | * sgs + | * chebyshev + | * iluk + | * ilut + | * icc + | * ict + | * amg + | * mgr + | * block + | * direct +amgNumSweeps integer 2 AMG smoother sweeps +amgSmootherType geosx_LinearSolverParameters_PreconditionerType gs | AMG smoother type. Valid options (not all may be supported by linear algebra package): + | * none + | * jacobi + | * gs + | * sgs + | * chebyshev + | * iluk + | * ilut + | * icc + | * ict + | * amg + | * mgr + | * block + | * direct +amgThreshold real64 0 AMG strength-of-connection threshold +directCheckResTol real64 1e-12 Tolerance used to check a direct solver solution +directColPerm geosx_LinearSolverParameters_Direct_ColPerm metis | How to permute the columns. Available options are: + | * none + | * MMD_AtplusA + | * MMD_AtA + | * colAMD + | * metis + | * parmetis +directEquil integer 1 Whether to scale the rows and columns of the matrix +directIterRef integer 1 Whether to perform iterative refinement +directParallel integer 1 Whether to use a parallel solver (instead of a serial one) +directReplTinyPivot integer 1 Whether to replace tiny pivots by sqrt(epsilon)*norm(A) +directRowPerm geosx_LinearSolverParameters_Direct_RowPerm mc64 | How to permute the rows. Available options are: + | * none + | * mc64 +iluFill integer 0 ILU(K) fill factor +iluThreshold real64 0 ILU(T) threshold factor +krylovAdaptiveTol integer 0 Use Eisenstat-Walker adaptive linear tolerance +krylovMaxIter integer 200 Maximum iterations allowed for an iterative solver +krylovMaxRestart integer 200 Maximum iterations before restart (GMRES only) +krylovTol real64 1e-06 | Relative convergence tolerance of the iterative method + | If the method converges, the iterative solution :math:`\mathsf{x}_k` is such that + | the relative residual norm satisfies: + | :math:`\left\lVert \mathsf{b} - \mathsf{A} \mathsf{x}_k \right\rVert_2` < ``krylovTol`` * :math:`\left\lVert\mathsf{b}\right\rVert_2` +krylovWeakestTol real64 0.001 Weakest-allowed tolerance for adaptive method +logLevel integer 0 Log level +preconditionerType geosx_LinearSolverParameters_PreconditionerType iluk | Preconditioner type. Available options are: + | * none + | * jacobi + | * gs + | * sgs + | * chebyshev + | * iluk + | * ilut + | * icc + | * ict + | * amg + | * mgr + | * block + | * direct +solverType geosx_LinearSolverParameters_SolverType direct | Linear solver type. Available options are: + | * direct + | * cg + | * gmres + | * fgmres + | * bicgstab + | * preconditioner +stopIfError integer 1 Whether to stop the simulation if the linear solver reports an error +=================== =============================================== ======= ======================================================================================================================================================================================================================================================================================================================= diff --git a/src/coreComponents/fileIO/schema/schema.xsd b/src/coreComponents/fileIO/schema/schema.xsd index 533c523576c..8687775fbe0 100644 --- a/src/coreComponents/fileIO/schema/schema.xsd +++ b/src/coreComponents/fileIO/schema/schema.xsd @@ -575,14 +575,38 @@ - - + + - - + + @@ -631,13 +655,15 @@ the relative residual norm satisfies: * jacobi * gs * sgs +* chebyshev * iluk * ilut * icc * ict * amg * mgr -* block--> +* block +* direct--> + + + + + @@ -660,11 +691,6 @@ the relative residual norm satisfies: - - - - - diff --git a/src/coreComponents/finiteVolume/CMakeLists.txt b/src/coreComponents/finiteVolume/CMakeLists.txt index f54e66bb021..8e25dae1081 100644 --- 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 index 136b04c732a..9a29d9da02b 100644 --- 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 index ad2a56a74d8..5aaed0f7790 100644 --- 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 index 8a7ed017950..21e18a0f481 100644 --- 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 index 64d2ed742b9..8cfb1b9b3a3 100644 --- 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 index 0be86180367..76b74dfe127 100644 --- 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 index c8df0bd85e0..7cf3bba9741 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.cpp @@ -25,6 +25,7 @@ #include "HypreUtils.hpp" #include +#include namespace geosx { @@ -447,53 +448,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, @@ -908,10 +862,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 index 3d633f2c939..2cdf944570d 100644 --- 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 index fafc68cc8ce..92c7ccc33ba 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.cpp @@ -74,6 +74,11 @@ HyprePreconditioner::HyprePreconditioner( LinearSolverParameters params, createILUT(); break; } + case LinearSolverParameters::PreconditionerType::direct: + { + createDirect(); + break; + } default: { GEOSX_ERROR( "Preconditioner type not supported in hypre interface: " << m_parameters.preconditionerType ); @@ -102,18 +107,23 @@ HYPRE_Int getHypreAMGCycleType( string const & type ) return typeMap.at( type ); } -HYPRE_Int getHypreAMGRelaxationType( string const & type ) +HYPRE_Int getHypreAMGRelaxationType( LinearSolverParameters::PreconditionerType const & type ) { - static std::map< string, HYPRE_Int > const typeMap = + static std::map< LinearSolverParameters::PreconditionerType, HYPRE_Int > const typeMap = { - { "jacobi", 0 }, + { LinearSolverParameters::PreconditionerType::jacobi, 0 }, + { LinearSolverParameters::PreconditionerType::gs, 6 }, + { LinearSolverParameters::PreconditionerType::chebyshev, 16 }, + + // TODO: re-enable if/when we add those options to PreconditionerType + // or find another way to pass L1/hybrid flags +#if 0 { "hybridForwardGaussSeidel", 3 }, { "hybridBackwardGaussSeidel", 4 }, { "hybridSymmetricGaussSeidel", 6 }, - { "gaussSeidel", 6 }, { "L1hybridSymmetricGaussSeidel", 8 }, - { "chebyshev", 16 }, { "L1jacobi", 18 }, +#endif }; GEOSX_LAI_ASSERT_MSG( typeMap.count( type ) > 0, "Unsupported Hypre AMG relaxation option: " << type ); @@ -140,20 +150,24 @@ void HyprePreconditioner::createAMG() // Set smoother to be used (other options available, see hypre's documentation) // (default "gaussSeidel", i.e. local symmetric Gauss-Seidel) - if( m_parameters.amg.smootherType.substr( 0, 3 ) == "ilu" ) + if( m_parameters.amg.smootherType == LinearSolverParameters::PreconditionerType::iluk ) { GEOSX_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetSmoothType( m_precond, 5 ) ); GEOSX_LAI_CHECK_ERROR( HYPRE_ILUSetType( m_precond, 0 ) ); + + // TODO: Re-enable disable ILU(1) option when we can properly pass level-of-fill for smoother +#if 0 if( m_parameters.amg.smootherType == "ilu1" ) { GEOSX_LAI_CHECK_ERROR( HYPRE_ILUSetLevelOfFill( m_precond, 1 ) ); } +#endif } else { GEOSX_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetRelaxType( m_precond, getHypreAMGRelaxationType( m_parameters.amg.smootherType ) ) ); - if( m_parameters.amg.smootherType == "chebyshev" ) + if( m_parameters.amg.smootherType == LinearSolverParameters::PreconditionerType::chebyshev ) { // Set order for Chebyshev smoother valid options 1, 2 (default), 3, 4) if( ( 0 < m_parameters.amg.numSweeps ) && ( m_parameters.amg.numSweeps < 5 ) ) @@ -166,7 +180,7 @@ void HyprePreconditioner::createAMG() // Set coarsest level solver // (by default for coarsest grid size above 5,000 superlu_dist is used) GEOSX_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetDSLUThreshold( m_precond, 5000 ) ); - if( m_parameters.amg.coarseType == "direct" ) + if( m_parameters.amg.coarseType == LinearSolverParameters::PreconditionerType::direct ) { GEOSX_LAI_CHECK_ERROR( hypre_BoomerAMGSetCycleRelaxType( m_precond, 9, 3 ) ); } @@ -542,6 +556,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 index f85b8a0443a..04ed9773482 100644 --- 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 index c2eab315c5b..c6c440d525a 100644 --- 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 index 591ce96b3ea..9a0e40a3344 100644 --- 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; diff --git a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.cpp b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.cpp index f267377826d..e3f1ca1a7fb 100644 --- a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.cpp @@ -314,7 +314,7 @@ void PetscMatrix::insert( globalIndex const rowIndex, void PetscMatrix::add( arraySlice1d< globalIndex const > const & rowIndices, arraySlice1d< globalIndex const > const & colIndices, - arraySlice2d< real64 const, 1 > const & values ) + arraySlice2d< real64 const, MatrixLayout::ROW_MAJOR > const & values ) { GEOSX_LAI_ASSERT( modifiable() ); GEOSX_LAI_CHECK_ERROR( MatSetValues( m_mat, @@ -328,7 +328,7 @@ void PetscMatrix::add( arraySlice1d< globalIndex const > const & rowIndices, void PetscMatrix::set( arraySlice1d< globalIndex const > const & rowIndices, arraySlice1d< globalIndex const > const & colIndices, - arraySlice2d< real64 const, 1 > const & values ) + arraySlice2d< real64 const, MatrixLayout::ROW_MAJOR > const & values ) { GEOSX_LAI_ASSERT( modifiable() ); GEOSX_LAI_CHECK_ERROR( MatSetValues( m_mat, @@ -342,7 +342,7 @@ void PetscMatrix::set( arraySlice1d< globalIndex const > const & rowIndices, void PetscMatrix::insert( arraySlice1d< globalIndex const > const & rowIndices, arraySlice1d< globalIndex const > const & colIndices, - arraySlice2d< real64 const, 1 > const & values ) + arraySlice2d< real64 const, MatrixLayout::ROW_MAJOR > const & values ) { GEOSX_LAI_ASSERT( insertable() ); GEOSX_LAI_CHECK_ERROR( MatSetValues( m_mat, @@ -354,54 +354,6 @@ void PetscMatrix::insert( arraySlice1d< globalIndex const > const & rowIndices, ADD_VALUES ) ); } -void PetscMatrix::add( arraySlice1d< globalIndex const > const & rowIndices, - arraySlice1d< globalIndex const > const & colIndices, - arraySlice2d< real64 const, 0 > const & values ) -{ - GEOSX_LAI_ASSERT( modifiable() ); - GEOSX_LAI_CHECK_ERROR( MatSetOption( m_mat, MAT_ROW_ORIENTED, PETSC_FALSE ) ); - GEOSX_LAI_CHECK_ERROR( MatSetValues( m_mat, - rowIndices.size(), - toPetscInt( rowIndices ), - colIndices.size(), - toPetscInt( colIndices ), - values.dataIfContiguous(), - ADD_VALUES ) ); - GEOSX_LAI_CHECK_ERROR( MatSetOption( m_mat, MAT_ROW_ORIENTED, PETSC_TRUE ) ); -} - -void PetscMatrix::set( arraySlice1d< globalIndex const > const & rowIndices, - arraySlice1d< globalIndex const > const & colIndices, - arraySlice2d< real64 const, 0 > const & values ) -{ - GEOSX_LAI_ASSERT( modifiable() ); - GEOSX_LAI_CHECK_ERROR( MatSetOption( m_mat, MAT_ROW_ORIENTED, PETSC_FALSE ) ); - GEOSX_LAI_CHECK_ERROR( MatSetValues( m_mat, - rowIndices.size(), - toPetscInt( rowIndices ), - colIndices.size(), - toPetscInt( colIndices ), - values.dataIfContiguous(), - INSERT_VALUES ) ); - GEOSX_LAI_CHECK_ERROR( MatSetOption( m_mat, MAT_ROW_ORIENTED, PETSC_TRUE ) ); -} - -void PetscMatrix::insert( arraySlice1d< globalIndex const > const & rowIndices, - arraySlice1d< globalIndex const > const & colIndices, - arraySlice2d< real64 const, 0 > const & values ) -{ - GEOSX_LAI_ASSERT( insertable() ); - GEOSX_LAI_CHECK_ERROR( MatSetOption( m_mat, MAT_ROW_ORIENTED, PETSC_FALSE ) ); - GEOSX_LAI_CHECK_ERROR( MatSetValues( m_mat, - rowIndices.size(), - toPetscInt( rowIndices ), - colIndices.size(), - toPetscInt( colIndices ), - values.dataIfContiguous(), - ADD_VALUES ) ); - GEOSX_LAI_CHECK_ERROR( MatSetOption( m_mat, MAT_ROW_ORIENTED, PETSC_TRUE ) ); -} - void PetscMatrix::add( globalIndex const * rowIndices, globalIndex const * colIndices, real64 const * values, @@ -661,19 +613,19 @@ localIndex PetscMatrix::maxRowLength() const { GEOSX_LAI_ASSERT( assembled() ); localIndex maxLocalLength = 0; - for( localIndex i = ilower(); i < iupper(); ++i ) + for( globalIndex i = ilower(); i < iupper(); ++i ) { maxLocalLength = std::max( maxLocalLength, globalRowLength( i ) ); } return MpiWrapper::Max( maxLocalLength, getComm() ); } -localIndex PetscMatrix::localRowLength( localIndex localRowIndex ) const +localIndex PetscMatrix::localRowLength( localIndex const localRowIndex ) const { return globalRowLength( getGlobalRowID( localRowIndex ) ); } -localIndex PetscMatrix::globalRowLength( globalIndex globalRowIndex ) const +localIndex PetscMatrix::globalRowLength( globalIndex const globalRowIndex ) const { GEOSX_LAI_ASSERT( assembled() ); PetscInt ncols; @@ -890,9 +842,9 @@ localIndex PetscMatrix::numLocalCols() const localIndex PetscMatrix::numLocalRows() const { GEOSX_LAI_ASSERT( created() ); - PetscInt low, high; - GEOSX_LAI_CHECK_ERROR( MatGetOwnershipRange( m_mat, &low, &high ) ); - return LvArray::integerConversion< localIndex >( high - low ); + PetscInt rows; + GEOSX_LAI_CHECK_ERROR( MatGetLocalSize( m_mat, &rows, nullptr ) ); + return LvArray::integerConversion< localIndex >( rows ); } MPI_Comm PetscMatrix::getComm() const diff --git a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.hpp b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.hpp index 54ca7fa66c6..a8fe1b3eb74 100644 --- a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.hpp @@ -166,18 +166,6 @@ class PetscMatrix final : public virtual LinearOperator< PetscVector >, 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/petsc/PetscPreconditioner.cpp b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscPreconditioner.cpp index e8ed988d354..c3ced37cf84 100644 --- a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscPreconditioner.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscPreconditioner.cpp @@ -18,8 +18,6 @@ #include "PetscPreconditioner.hpp" -#include "linearAlgebra/utilities/LAIHelperFunctions.hpp" - #include namespace geosx @@ -37,7 +35,7 @@ PetscPreconditioner::~PetscPreconditioner() clear(); } -void CreatePetscAMG( LinearSolverParameters const & params, PC precond ) +void CreatePetscAMG( LinearSolverParameters const & params, PC const precond ) { // Default options only for the moment GEOSX_LAI_CHECK_ERROR( PCSetType( precond, PCGAMG ) ); @@ -45,32 +43,32 @@ void CreatePetscAMG( LinearSolverParameters const & params, PC precond ) // TODO: need someone familiar with PETSc to take a look at this #if 0 - GEOSX_LAI_CHECK_ERROR( PCSetType( prec, PCHMG ) ); - GEOSX_LAI_CHECK_ERROR( PCHMGSetInnerPCType( prec, PCGAMG ) ); + GEOSX_LAI_CHECK_ERROR( PCSetType( precond, PCHMG ) ); + GEOSX_LAI_CHECK_ERROR( PCHMGSetInnerPCType( precond, PCGAMG ) ); // Set maximum number of multigrid levels - if( m_parameters.amg.maxLevels > 0 ) + if( params.amg.maxLevels > 0 ) { - GEOSX_LAI_CHECK_ERROR( PCMGSetLevels( prec, - LvArray::integerConversion< PetscInt >( m_parameters.amg.maxLevels ), + GEOSX_LAI_CHECK_ERROR( PCMGSetLevels( precond, + LvArray::integerConversion< PetscInt >( params.amg.maxLevels ), nullptr ) ); } // Set the number of sweeps - if( m_parameters.amg.numSweeps > 1 ) + if( params.amg.numSweeps > 1 ) { - GEOSX_LAI_CHECK_ERROR( PCMGSetNumberSmooth( prec, - LvArray::integerConversion< PetscInt >( m_parameters.amg.numSweeps ) ) ); + GEOSX_LAI_CHECK_ERROR( PCMGSetNumberSmooth( precond, + LvArray::integerConversion< PetscInt >( params.amg.numSweeps ) ) ); } // Set type of cycle (1: V-cycle (default); 2: W-cycle) - if( m_parameters.amg.cycleType == "V" ) + if( params.amg.cycleType == "V" ) { - GEOSX_LAI_CHECK_ERROR( PCMGSetCycleType( prec, PC_MG_CYCLE_V ) ); + GEOSX_LAI_CHECK_ERROR( PCMGSetCycleType( precond, PC_MG_CYCLE_V ) ); } - else if( m_parameters.amg.cycleType == "W" ) + else if( params.amg.cycleType == "W" ) { - GEOSX_LAI_CHECK_ERROR( PCMGSetCycleType( prec, PC_MG_CYCLE_W ) ); + GEOSX_LAI_CHECK_ERROR( PCMGSetCycleType( precond, PC_MG_CYCLE_W ) ); } // Set smoother to be used (for all levels) @@ -78,55 +76,67 @@ void CreatePetscAMG( LinearSolverParameters const & params, PC precond ) PetscInt l; KSP smoother; PC smootherPC; - GEOSX_LAI_CHECK_ERROR( PCMGGetLevels( prec, &numLevels ) ); + GEOSX_LAI_CHECK_ERROR( PCMGGetLevels( precond, &numLevels ) ); GEOSX_LOG_RANK_VAR( numLevels ); for( l = 0; l < numLevels; ++l ) { - GEOSX_LAI_CHECK_ERROR( PCMGGetSmoother( prec, l, &smoother ) ); + GEOSX_LAI_CHECK_ERROR( PCMGGetSmoother( precond, l, &smoother ) ); GEOSX_LAI_CHECK_ERROR( KSPSetType( smoother, KSPRICHARDSON ) ); GEOSX_LAI_CHECK_ERROR( KSPGetPC( smoother, &smootherPC ) ); - if( m_parameters.amg.smootherType == "jacobi" ) - { - GEOSX_LAI_CHECK_ERROR( PCSetType( smootherPC, PCJACOBI ) ); - } - else if( m_parameters.amg.smootherType == "gaussSeidel" ) + switch( params.amg.smootherType ) { - GEOSX_LAI_CHECK_ERROR( PCSetType( smootherPC, PCSOR ) ); - } - else if( m_parameters.amg.smootherType.substr( 0, 3 ) == "ilu" ) - { - // Set up additive Schwartz preconditioner - GEOSX_LAI_CHECK_ERROR( PCSetType( smootherPC, PCASM ) ); - GEOSX_LAI_CHECK_ERROR( PCASMSetOverlap( smootherPC, 0 ) ); - GEOSX_LAI_CHECK_ERROR( PCASMSetType( smootherPC, PC_ASM_RESTRICT ) ); - // GEOSX_LAI_CHECK_ERROR( PCSetUp( smootherPC ) ); - - // Get local preconditioning context - KSP * ksp_local; - PetscInt n_local, first_local; - GEOSX_LAI_CHECK_ERROR( PCASMGetSubKSP( smootherPC, &n_local, &first_local, &ksp_local ) ); - - // Sanity checks - GEOSX_LAI_ASSERT_EQ( n_local, 1 ); - GEOSX_LAI_ASSERT_EQ( first_local, MpiWrapper::Comm_rank( MPI_COMM_GEOSX ) ); - - // Set up local block ILU preconditioner - PC prec_local; - GEOSX_LAI_CHECK_ERROR( KSPSetType( ksp_local[0], KSPPREONLY ) ); - GEOSX_LAI_CHECK_ERROR( KSPGetPC( ksp_local[0], &prec_local ) ); - GEOSX_LAI_CHECK_ERROR( PCSetType( prec_local, PCILU ) ); - if( m_parameters.amg.smootherType == "ilu1" ) + case LinearSolverParameters::PreconditionerType::jacobi: { - GEOSX_LAI_CHECK_ERROR( PCFactorSetLevels( prec_local, 1 ) ); + GEOSX_LAI_CHECK_ERROR( PCSetType( smootherPC, PCJACOBI ) ); + break; } - else + case LinearSolverParameters::PreconditionerType::gs: + { + GEOSX_LAI_CHECK_ERROR( PCSetType( smootherPC, PCSOR ) ); + break; + } + case LinearSolverParameters::PreconditionerType::iluk: { - GEOSX_LAI_CHECK_ERROR( PCFactorSetLevels( prec_local, 0 ) ); + // Set up additive Schwartz preconditioner + GEOSX_LAI_CHECK_ERROR( PCSetType( smootherPC, PCASM ) ); + GEOSX_LAI_CHECK_ERROR( PCASMSetOverlap( smootherPC, 0 ) ); + GEOSX_LAI_CHECK_ERROR( PCASMSetType( smootherPC, PC_ASM_RESTRICT ) ); + // GEOSX_LAI_CHECK_ERROR( PCSetUp( smootherPC ) ); + + // Get local preconditioning context + KSP * ksp_local; + PetscInt n_local, first_local; + GEOSX_LAI_CHECK_ERROR( PCASMGetSubKSP( smootherPC, &n_local, &first_local, &ksp_local ) ); + + // Sanity checks + GEOSX_LAI_ASSERT_EQ( n_local, 1 ); + GEOSX_LAI_ASSERT_EQ( first_local, MpiWrapper::Comm_rank( MPI_COMM_GEOSX ) ); + + // Set up local block ILU preconditioner + PC prec_local; + GEOSX_LAI_CHECK_ERROR( KSPSetType( ksp_local[0], KSPPREONLY ) ); + GEOSX_LAI_CHECK_ERROR( KSPGetPC( ksp_local[0], &prec_local ) ); + GEOSX_LAI_CHECK_ERROR( PCSetType( prec_local, PCILU ) ); + + // TODO: re-enable when we can properly pass level-of-fill in smoother +#if 0 + if( params.amg.smootherType == "ilu1" ) + { + GEOSX_LAI_CHECK_ERROR( PCFactorSetLevels( prec_local, 1 ) ); + } + else +#endif + { + GEOSX_LAI_CHECK_ERROR( PCFactorSetLevels( prec_local, 0 ) ); + } + // GEOSX_LAI_CHECK_ERROR( PCSetUpOnBlocks( smootherPC ) ); + break; } - // GEOSX_LAI_CHECK_ERROR( PCSetUpOnBlocks( smootherPC ) ); + default: + GEOSX_ERROR( "Smoother type not supported in PETSc/AMG: " << params.amg.smootherType ); } } @@ -161,7 +171,7 @@ PCType getPetscSmootherType( LinearSolverParameters::PreconditionerType const & return typeMap.at( type ); } -void CreatePetscSmoother( LinearSolverParameters const & params, PC precond ) +void CreatePetscSmoother( LinearSolverParameters const & params, PC const precond ) { // Set up additive Schwartz outer preconditioner GEOSX_LAI_CHECK_ERROR( PCSetType( precond, PCASM ) ); @@ -183,6 +193,13 @@ void CreatePetscSmoother( LinearSolverParameters const & params, PC precond ) GEOSX_LAI_CHECK_ERROR( PCFactorSetLevels( prec_local, params.ilu.fill ) ); } +void CreatePetscDirect( LinearSolverParameters const & GEOSX_UNUSED_PARAM( params ), PC const precond ) +{ + GEOSX_LAI_CHECK_ERROR( PCSetType( precond, PCLU ) ); + GEOSX_LAI_CHECK_ERROR( PCFactorSetMatSolverType( precond, MATSOLVERSUPERLU_DIST ) ); + GEOSX_LAI_CHECK_ERROR( PCSetUp( precond ) ); +} + void PetscPreconditioner::compute( PetscMatrix const & mat ) { Base::compute( mat ); @@ -222,6 +239,11 @@ void PetscPreconditioner::compute( PetscMatrix const & mat ) CreatePetscSmoother( m_parameters, m_precond ); break; } + case LinearSolverParameters::PreconditionerType::direct: + { + CreatePetscDirect( m_parameters, m_precond ); + break; + } default: { GEOSX_ERROR( "Preconditioner type not supported in PETSc interface: " << m_parameters.preconditionerType ); diff --git a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscVector.cpp b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscVector.cpp index 94eaf698118..3861d8aaa5f 100644 --- a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscVector.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscVector.cpp @@ -107,7 +107,7 @@ void PetscVector::createWithGlobalSize( globalIndex const globalSize, MPI_Comm c GEOSX_LAI_CHECK_ERROR( VecSetSizes( m_vec, PETSC_DECIDE, globalSize ) ); } -void PetscVector::create( arrayView1d< real64 const > const & localValues, MPI_Comm const & comm ) +void PetscVector::createWithLocalValues( arrayView1d< real64 > const & localValues, MPI_Comm const & comm ) { GEOSX_LAI_ASSERT( closed() ); reset(); @@ -447,7 +447,7 @@ globalIndex PetscVector::globalSize() const GEOSX_LAI_ASSERT( created() ); PetscInt size; GEOSX_LAI_CHECK_ERROR( VecGetSize( m_vec, &size ) ); - return size; + return LvArray::integerConversion< globalIndex >( size ); } localIndex PetscVector::localSize() const @@ -455,7 +455,7 @@ localIndex PetscVector::localSize() const GEOSX_LAI_ASSERT( created() ); PetscInt size; GEOSX_LAI_CHECK_ERROR( VecGetLocalSize( m_vec, &size ) ); - return size; + return LvArray::integerConversion< localIndex >( size ); } localIndex PetscVector::getLocalRowID( globalIndex const globalRow ) const diff --git a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscVector.hpp b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscVector.hpp index 044f04c13fa..abfbc403c78 100644 --- a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscVector.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscVector.hpp @@ -108,8 +108,8 @@ class PetscVector final : private VectorBase< PetscVector > 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; diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.cpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.cpp index 6738dcd0953..7b6c6f6403b 100644 --- a/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.cpp @@ -30,13 +30,6 @@ #include #include -#ifdef GEOSX_USE_MPI -#include -#else -#include -using Epetra_MpiComm = Epetra_SerialComm; -#endif - namespace geosx { @@ -71,10 +64,10 @@ void EpetraMatrix::createWithGlobalSize( globalIndex const globalRows, m_dst_map = std::make_unique< Epetra_Map >( globalRows, 0, - Epetra_MpiComm( MPI_PARAM( comm ) ) ); + EpetraComm( MPI_PARAM( comm ) ) ); m_src_map = std::make_unique< Epetra_Map >( globalCols, 0, - Epetra_MpiComm( MPI_PARAM( comm ) ) ); + EpetraComm( MPI_PARAM( comm ) ) ); m_matrix = std::make_unique< Epetra_FECrsMatrix >( Copy, *m_dst_map, LvArray::integerConversion< int >( maxEntriesPerRow ), @@ -96,11 +89,11 @@ void EpetraMatrix::createWithLocalSize( localIndex const localRows, m_dst_map = std::make_unique< Epetra_Map >( LvArray::integerConversion< globalIndex >( -1 ), LvArray::integerConversion< int >( localRows ), 0, - Epetra_MpiComm( MPI_PARAM( comm ) ) ); + EpetraComm( MPI_PARAM( comm ) ) ); m_src_map = std::make_unique< Epetra_Map >( LvArray::integerConversion< globalIndex >( -1 ), LvArray::integerConversion< int >( localCols ), 0, - Epetra_MpiComm( MPI_PARAM( comm ) ) ); + EpetraComm( MPI_PARAM( comm ) ) ); m_matrix = std::make_unique< Epetra_FECrsMatrix >( Copy, *m_dst_map, LvArray::integerConversion< int >( maxEntriesPerRow ), @@ -277,45 +270,6 @@ void EpetraMatrix::insert( arraySlice1d< globalIndex const > const & rowIndices, Epetra_FECrsMatrix::ROW_MAJOR ) ); } -void EpetraMatrix::add( arraySlice1d< globalIndex const > const & rowIndices, - arraySlice1d< globalIndex const > const & colIndices, - arraySlice2d< real64 const, MatrixLayout::COL_MAJOR > const & values ) -{ - GEOSX_LAI_ASSERT( modifiable() ); - GEOSX_LAI_CHECK_ERROR( m_matrix->SumIntoGlobalValues( LvArray::integerConversion< int >( rowIndices.size() ), - toEpetraLongLong( rowIndices ), - LvArray::integerConversion< int >( colIndices.size() ), - toEpetraLongLong( colIndices ), - values.dataIfContiguous(), - Epetra_FECrsMatrix::COLUMN_MAJOR ) ); -} - -void EpetraMatrix::set( arraySlice1d< globalIndex const > const & rowIndices, - arraySlice1d< globalIndex const > const & colIndices, - arraySlice2d< real64 const, MatrixLayout::COL_MAJOR > const & values ) -{ - GEOSX_LAI_ASSERT( modifiable() ); - GEOSX_LAI_CHECK_ERROR( m_matrix->ReplaceGlobalValues( LvArray::integerConversion< int >( rowIndices.size() ), - toEpetraLongLong( rowIndices ), - LvArray::integerConversion< int >( colIndices.size() ), - toEpetraLongLong( colIndices ), - values.dataIfContiguous(), - Epetra_FECrsMatrix::COLUMN_MAJOR ) ); -} - -void EpetraMatrix::insert( arraySlice1d< globalIndex const > const & rowIndices, - arraySlice1d< globalIndex const > const & colIndices, - arraySlice2d< real64 const, MatrixLayout::COL_MAJOR > const & values ) -{ - GEOSX_LAI_ASSERT( insertable() ); - GEOSX_LAI_CHECK_ERROR_NNEG( m_matrix->InsertGlobalValues( LvArray::integerConversion< int >( rowIndices.size() ), - toEpetraLongLong( rowIndices ), - LvArray::integerConversion< int >( colIndices.size() ), - toEpetraLongLong( colIndices ), - values.dataIfContiguous(), - Epetra_FECrsMatrix::COLUMN_MAJOR ) ); -} - void EpetraMatrix::add( globalIndex const * rowIndices, globalIndex const * colIndices, real64 const * values, @@ -594,7 +548,7 @@ void EpetraMatrix::addDiagonal( EpetraVector const & src ) localIndex EpetraMatrix::maxRowLength() const { GEOSX_LAI_ASSERT( assembled() ); - return m_matrix->MaxNumEntries(); + return m_matrix->GlobalMaxNumEntries(); } localIndex EpetraMatrix::localRowLength( localIndex localRowIndex ) const @@ -782,13 +736,13 @@ globalIndex EpetraMatrix::getGlobalRowID( localIndex const index ) const localIndex EpetraMatrix::numLocalCols() const { GEOSX_LAI_ASSERT( created() ); - return m_src_map->NumMyElements(); + return LvArray::integerConversion< localIndex >( m_src_map->NumMyElements() ); } localIndex EpetraMatrix::numLocalRows() const { GEOSX_LAI_ASSERT( created() ); - return m_matrix->RowMap().NumMyElements(); + return LvArray::integerConversion< localIndex >( m_dst_map->NumMyElements() ); } MPI_Comm EpetraMatrix::getComm() const diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.hpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.hpp index bd3e2b812ef..c2be798fe6d 100644 --- a/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.hpp @@ -19,10 +19,9 @@ #ifndef GEOSX_LINEARALGEBRA_INTERFACES_EPETRAMATRIX_HPP_ #define GEOSX_LINEARALGEBRA_INTERFACES_EPETRAMATRIX_HPP_ -#include "common/DataTypes.hpp" -#include "linearAlgebra/interfaces/trilinos/EpetraVector.hpp" #include "linearAlgebra/interfaces/LinearOperator.hpp" #include "linearAlgebra/interfaces/MatrixBase.hpp" +#include "linearAlgebra/interfaces/trilinos/EpetraVector.hpp" class Epetra_Map; @@ -157,18 +156,6 @@ class EpetraMatrix final : public virtual LinearOperator< EpetraVector >, 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/trilinos/EpetraUtils.hpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraUtils.hpp index dfc1ab80a2e..c12da401bed 100644 --- a/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraUtils.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraUtils.hpp @@ -20,6 +20,18 @@ #include "common/DataTypes.hpp" +#ifdef GEOSX_USE_MPI +#include + +/// Alias for specific EpetraComm implementation used. +using EpetraComm = Epetra_MpiComm; +#else +#include + +/// Alias for specific EpetraComm implementation used. +using EpetraComm = Epetra_SerialComm; +#endif + namespace geosx { diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraVector.cpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraVector.cpp index 95cb6387cfe..cad81e343b3 100644 --- a/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraVector.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraVector.cpp @@ -25,13 +25,6 @@ #include #include -#ifdef GEOSX_USE_MPI -#include -#else -#include -typedef Epetra_SerialComm Epetra_MpiComm; -#endif - namespace geosx { @@ -101,7 +94,7 @@ void EpetraVector::createWithLocalSize( localIndex const localSize, Epetra_Map const map( LvArray::integerConversion< long long >( -1 ), LvArray::integerConversion< int >( localSize ), 0, - Epetra_MpiComm( MPI_PARAM( comm ) ) ); + EpetraComm( MPI_PARAM( comm ) ) ); m_vector = std::make_unique< Epetra_FEVector >( map ); } @@ -112,12 +105,12 @@ void EpetraVector::createWithGlobalSize( globalIndex const globalSize, GEOSX_LAI_ASSERT_GE( globalSize, 0 ); Epetra_Map const map( LvArray::integerConversion< long long >( globalSize ), 0, - Epetra_MpiComm( MPI_PARAM( comm ) ) ); + EpetraComm( MPI_PARAM( comm ) ) ); m_vector = std::make_unique< Epetra_FEVector >( map ); } -void EpetraVector::create( arrayView1d< real64 const > const & localValues, - MPI_Comm const & MPI_PARAM( comm ) ) +void EpetraVector::createWithLocalValues( arrayView1d< real64 > const & localValues, + MPI_Comm const & MPI_PARAM( comm ) ) { GEOSX_LAI_ASSERT( closed() ); @@ -127,10 +120,10 @@ void EpetraVector::create( arrayView1d< real64 const > const & localValues, Epetra_Map const map( LvArray::integerConversion< long long >( -1 ), localSize, 0, - Epetra_MpiComm( MPI_PARAM( comm ) ) ); - m_vector = std::make_unique< Epetra_FEVector >( Copy, + EpetraComm( MPI_PARAM( comm ) ) ); + m_vector = std::make_unique< Epetra_FEVector >( View, map, - const_cast< real64 * >( localValues.data() ), + localValues.data(), localSize, 1 ); } diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraVector.hpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraVector.hpp index dc15a85463b..edb06a21ff5 100644 --- a/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraVector.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraVector.hpp @@ -96,8 +96,8 @@ class EpetraVector final : private VectorBase< EpetraVector > 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; diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/TpetraMatrix.cpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TpetraMatrix.cpp new file mode 100644 index 00000000000..324ea93f14f --- /dev/null +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TpetraMatrix.cpp @@ -0,0 +1,795 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2019 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2019 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2019 Total, S.A + * Copyright (c) 2019- GEOSX Contributors + * All right reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file TpetraMatrix.cpp + */ + +#include "TpetraMatrix.hpp" + +#include "codingUtilities/Utilities.hpp" +#include "linearAlgebra/interfaces/trilinos/TpetraUtils.hpp" + +#include +#include +#include +#include +#include + +namespace geosx +{ + +TpetraMatrix::TpetraMatrix() + : LinearOperator(), + MatrixBase() +{ } + +TpetraMatrix::TpetraMatrix( TpetraMatrix const & src ) + : TpetraMatrix() +{ + GEOSX_LAI_ASSERT( src.ready() ); + m_matrix = std::make_unique< Tpetra_CrsMatrix >( *src.m_matrix ); + m_src_map = std::make_unique< Tpetra_Map >( *m_matrix->getDomainMap() ); + m_dst_map = std::make_unique< Tpetra_Map >( *m_matrix->getRangeMap() ); + m_assembled = true; +} + +TpetraMatrix::~TpetraMatrix() = default; + +void TpetraMatrix::createWithGlobalSize( globalIndex const globalRows, + globalIndex const globalCols, + localIndex const maxEntriesPerRow, + MPI_Comm const & comm ) +{ + GEOSX_LAI_ASSERT( closed() ); + GEOSX_LAI_ASSERT_GE( globalRows, 0 ); + GEOSX_LAI_ASSERT_GE( globalCols, 0 ); + GEOSX_LAI_ASSERT_GE( maxEntriesPerRow, 0 ); + + reset(); + + Teuchos::RCP< Tpetra_Comm > tpetraComm( new Tpetra_Comm( MPI_PARAM( comm ) ) ); + m_dst_map = std::make_unique< Tpetra_Map >( LvArray::integerConversion< Tpetra::global_size_t >( globalRows ), + 0, + tpetraComm ); + m_src_map = std::make_unique< Tpetra_Map >( LvArray::integerConversion< Tpetra::global_size_t >( globalCols ), + 0, + tpetraComm ); + m_matrix = std::make_unique< Tpetra_CrsMatrix >( Teuchos::RCP< Tpetra_Map >( m_dst_map.get(), false ), + LvArray::integerConversion< size_t >( maxEntriesPerRow ) ); +} + +void TpetraMatrix::createWithLocalSize( localIndex const localRows, + localIndex const localCols, + localIndex const maxEntriesPerRow, + MPI_Comm const & comm ) +{ + GEOSX_LAI_ASSERT( closed() ); + GEOSX_LAI_ASSERT_GE( localRows, 0 ); + GEOSX_LAI_ASSERT_GE( localCols, 0 ); + GEOSX_LAI_ASSERT_GE( maxEntriesPerRow, 0 ); + + reset(); + + Teuchos::RCP< Tpetra_Comm > tpetraComm( new Tpetra_Comm( MPI_PARAM( comm ) ) ); + m_dst_map = std::make_unique< Tpetra_Map >( Teuchos::OrdinalTraits< Tpetra::global_size_t >::invalid(), + LvArray::integerConversion< size_t >( localRows ), + 0, + tpetraComm ); + m_src_map = std::make_unique< Tpetra_Map >( Teuchos::OrdinalTraits< Tpetra::global_size_t >::invalid(), + LvArray::integerConversion< size_t >( localCols ), + 0, + tpetraComm ); + m_matrix = std::make_unique< Tpetra_CrsMatrix >( Teuchos::RCP< Tpetra_Map >( m_dst_map.get(), false ), + LvArray::integerConversion< size_t >( maxEntriesPerRow ) ); +} + +bool TpetraMatrix::created() const +{ + return bool(m_matrix); +} + +void TpetraMatrix::reset() +{ + MatrixBase::reset(); + m_matrix.reset(); + m_dst_map.reset(); + m_src_map.reset(); +} + +void TpetraMatrix::set( real64 const value ) +{ + GEOSX_LAI_ASSERT( ready() ); + open(); + m_matrix->setAllToScalar( value ); + close(); +} + +void TpetraMatrix::zero() +{ + set( 0 ); +} + +void TpetraMatrix::open() +{ + GEOSX_LAI_ASSERT( created() && closed() ); + m_matrix->resumeFill(); + m_closed = false; +} + +void TpetraMatrix::close() +{ + GEOSX_LAI_ASSERT( !closed() ); + // TODO: pass a Teuchos::ParameterList( { "No Nonlocal Changes", true } ) + m_matrix->fillComplete( Teuchos::RCP< Tpetra_Map >( m_src_map.get(), false ), + Teuchos::RCP< Tpetra_Map >( m_dst_map.get(), false ) ); + m_assembled = true; + m_closed = true; +} + +void TpetraMatrix::add( globalIndex const rowIndex, + globalIndex const colIndex, + real64 const value ) +{ + GEOSX_LAI_ASSERT( modifiable() ); + localIndex const numEntries = m_matrix->sumIntoGlobalValues( rowIndex, 1, &value, &colIndex ); + GEOSX_LAI_ASSERT_EQ( numEntries, 1 ); +} + +void TpetraMatrix::set( globalIndex const rowIndex, + globalIndex const colIndex, + real64 const value ) +{ + GEOSX_LAI_ASSERT( modifiable() ); + localIndex const numEntries = m_matrix->replaceGlobalValues( rowIndex, 1, &value, &colIndex ); + GEOSX_LAI_ASSERT_EQ( numEntries, 1 ); +} + +void TpetraMatrix::insert( globalIndex const rowIndex, + globalIndex const colIndex, + real64 const value ) +{ + GEOSX_LAI_ASSERT( insertable() ); + m_matrix->insertGlobalValues( rowIndex, 1, &value, &colIndex ); +} + +void TpetraMatrix::add( globalIndex const rowIndex, + globalIndex const * colIndices, + real64 const * values, + localIndex const size ) +{ + GEOSX_LAI_ASSERT( modifiable() ); + localIndex const numEntries = m_matrix->sumIntoGlobalValues( rowIndex, size, values, colIndices ); + GEOSX_LAI_ASSERT_EQ( numEntries, size ); +} + +void TpetraMatrix::set( globalIndex const rowIndex, + globalIndex const * colIndices, + real64 const * values, + localIndex const size ) +{ + GEOSX_LAI_ASSERT( modifiable() ); + localIndex const numEntries = m_matrix->replaceGlobalValues( rowIndex, size, values, colIndices ); + GEOSX_LAI_ASSERT_EQ( numEntries, size ); +} + +void TpetraMatrix::insert( globalIndex const rowIndex, + globalIndex const * colIndices, + real64 const * values, + localIndex const size ) +{ + GEOSX_LAI_ASSERT( insertable() ); + m_matrix->insertGlobalValues( rowIndex, size, values, colIndices ); +} + +void TpetraMatrix::add( globalIndex const rowIndex, + arraySlice1d< globalIndex const > const & colIndices, + arraySlice1d< real64 const > const & values ) +{ + GEOSX_LAI_ASSERT( modifiable() ); + localIndex const numEntries = m_matrix->sumIntoGlobalValues( rowIndex, colIndices.size(), values, colIndices ); + GEOSX_LAI_ASSERT_EQ( numEntries, colIndices.size() ); +} + +void TpetraMatrix::set( globalIndex const rowIndex, + arraySlice1d< globalIndex const > const & colIndices, + arraySlice1d< real64 const > const & values ) +{ + GEOSX_LAI_ASSERT( modifiable() ); + localIndex const numEntries = m_matrix->replaceGlobalValues( rowIndex, colIndices.size(), values, colIndices ); + GEOSX_LAI_ASSERT_EQ( numEntries, colIndices.size() ); +} + +void TpetraMatrix::insert( globalIndex const rowIndex, + arraySlice1d< globalIndex const > const & colIndices, + arraySlice1d< real64 const > const & values ) +{ + GEOSX_LAI_ASSERT( insertable() ); + m_matrix->insertGlobalValues( rowIndex, colIndices.size(), values, colIndices ); +} + +void TpetraMatrix::add( arraySlice1d< globalIndex const > const & rowIndices, + arraySlice1d< globalIndex const > const & colIndices, + arraySlice2d< real64 const, MatrixLayout::ROW_MAJOR > const & values ) +{ + GEOSX_LAI_ASSERT( modifiable() ); + for( localIndex i = 0; i < rowIndices.size(); ++i ) + { + add( rowIndices[i], colIndices, values[i] ); + } +} + +void TpetraMatrix::set( arraySlice1d< globalIndex const > const & rowIndices, + arraySlice1d< globalIndex const > const & colIndices, + arraySlice2d< real64 const, MatrixLayout::ROW_MAJOR > const & values ) +{ + GEOSX_LAI_ASSERT( modifiable() ); + for( localIndex i = 0; i < rowIndices.size(); ++i ) + { + set( rowIndices[i], colIndices, values[i] ); + } +} + +void TpetraMatrix::insert( arraySlice1d< globalIndex const > const & rowIndices, + arraySlice1d< globalIndex const > const & colIndices, + arraySlice2d< real64 const, MatrixLayout::ROW_MAJOR > const & values ) +{ + for( localIndex i = 0; i < rowIndices.size(); ++i ) + { + insert( rowIndices[i], colIndices, values[i] ); + } +} + +void TpetraMatrix::add( globalIndex const * rowIndices, + globalIndex const * colIndices, + real64 const * values, + localIndex const numRows, + localIndex const numCols ) +{ + for( localIndex i = 0; i < numRows; ++i ) + { + add( rowIndices[i], colIndices, values + i * numCols, numCols ); + } +} + +void TpetraMatrix::set( globalIndex const * rowIndices, + globalIndex const * colIndices, + real64 const * values, + localIndex const numRows, + localIndex const numCols ) +{ + for( localIndex i = 0; i < numRows; ++i ) + { + set( rowIndices[i], colIndices, values + i * numCols, numCols ); + } +} + +void TpetraMatrix::insert( globalIndex const * rowIndices, + globalIndex const * colIndices, + real64 const * values, + localIndex const numRows, + localIndex const numCols ) +{ + for( localIndex i = 0; i < numRows; ++i ) + { + insert( rowIndices[i], colIndices, values + i * numCols, numCols ); + } +} + +void TpetraMatrix::apply( TpetraVector const & src, + TpetraVector & dst ) const +{ + GEOSX_LAI_ASSERT( ready() ); + GEOSX_LAI_ASSERT( src.ready() ); + GEOSX_LAI_ASSERT( dst.ready() ); + GEOSX_LAI_ASSERT_EQ( numGlobalRows(), dst.globalSize() ); + GEOSX_LAI_ASSERT_EQ( numGlobalCols(), src.globalSize() ); + + m_matrix->apply( src.unwrapped(), + dst.unwrapped(), + Teuchos::ETransp::NO_TRANS, + 1.0, + 0.0 ); +} + +void TpetraMatrix::applyTranspose( TpetraVector const & src, + TpetraVector & dst ) const +{ + GEOSX_LAI_ASSERT( ready() ); + GEOSX_LAI_ASSERT( src.ready() ); + GEOSX_LAI_ASSERT( dst.ready() ); + GEOSX_LAI_ASSERT_EQ( numGlobalCols(), dst.globalSize() ); + GEOSX_LAI_ASSERT_EQ( numGlobalRows(), src.globalSize() ); + + m_matrix->apply( src.unwrapped(), + dst.unwrapped(), + Teuchos::ETransp::TRANS, + 1.0, + 0.0 ); +} + +void TpetraMatrix::multiply( TpetraMatrix const & src, + TpetraMatrix & dst ) const +{ + this->multiply( false, src, false, dst ); +} + +void TpetraMatrix::leftMultiplyTranspose( TpetraMatrix const & src, + TpetraMatrix & dst ) const +{ + this->multiply( true, src, false, dst ); +} + +void TpetraMatrix::rightMultiplyTranspose( TpetraMatrix const & src, + TpetraMatrix & dst ) const +{ + src.multiply( false, *this, true, dst ); +} + +void TpetraMatrix::create( Tpetra_CrsMatrix const & src ) +{ + GEOSX_LAI_ASSERT( closed() ); + reset(); + + m_matrix = std::make_unique< Tpetra_CrsMatrix >( src ); + m_src_map = std::make_unique< Tpetra_Map >( *m_matrix->getDomainMap() ); + m_dst_map = std::make_unique< Tpetra_Map >( *m_matrix->getRangeMap() ); + m_assembled = true; +} + +void TpetraMatrix::multiplyRAP( TpetraMatrix const & R, + TpetraMatrix const & P, + TpetraMatrix & dst ) const +{ + // TODO: Tpetra::TripleMatrixMultiply::MultiplyRAP gives wrong results, investigate +#if 0 + GEOSX_LAI_ASSERT( ready() ); + GEOSX_LAI_ASSERT( R.ready() ); + GEOSX_LAI_ASSERT( P.ready() ); + GEOSX_LAI_ASSERT_EQ( R.numLocalCols(), numLocalRows() ); + GEOSX_LAI_ASSERT_EQ( P.numLocalRows(), numLocalCols() ); + + // TODO: estimate num nonzeros per row? or it does not matter since matrix is recreated anyway? + dst.createWithLocalSize( R.numLocalRows(), P.numLocalCols(), 1, getComm() ); + + Tpetra::TripleMatrixMultiply::MultiplyRAP( R.unwrapped(), + false, + *m_matrix, + false, + P.unwrapped(), + false, + dst.unwrapped() ); + dst.m_assembled = true; +#else + MatrixBase::multiplyRAP( R, P, dst ); +#endif +} + +void TpetraMatrix::multiplyPtAP( TpetraMatrix const & P, + TpetraMatrix & dst ) const +{ + GEOSX_LAI_ASSERT( ready() ); + GEOSX_LAI_ASSERT( P.ready() ); + GEOSX_LAI_ASSERT_EQ( P.numLocalRows(), numLocalRows() ); + GEOSX_LAI_ASSERT_EQ( P.numLocalRows(), numLocalCols() ); + + // TODO: estimate num nonzeros per row? or it does not matter since matrix is recreated anyway? + dst.createWithLocalSize( P.numLocalCols(), P.numLocalCols(), 1, getComm() ); + + Tpetra::TripleMatrixMultiply::MultiplyRAP( P.unwrapped(), + true, + *m_matrix, + false, + P.unwrapped(), + false, + dst.unwrapped() ); + dst.m_assembled = true; +} + +void TpetraMatrix::gemv( real64 const alpha, + TpetraVector const & x, + real64 const beta, + TpetraVector & y, + bool useTranspose ) const +{ + GEOSX_LAI_ASSERT( ready() ); + m_matrix->apply( x.unwrapped(), + y.unwrapped(), + useTranspose ? Teuchos::ETransp::TRANS : Teuchos::ETransp::NO_TRANS, + alpha, + beta ); +} + +void TpetraMatrix::scale( real64 const scalingFactor ) +{ + GEOSX_LAI_ASSERT( ready() ); + + if( isEqual( scalingFactor, 1.0 ) ) + { + return; + } + + open(); + m_matrix->scale( scalingFactor ); + close(); +} + +void TpetraMatrix::leftScale( TpetraVector const & vec ) +{ + GEOSX_LAI_ASSERT( ready() ); + m_matrix->leftScale( vec.unwrapped() ); +} + +void TpetraMatrix::rightScale( TpetraVector const & vec ) +{ + GEOSX_LAI_ASSERT( ready() ); + m_matrix->rightScale( vec.unwrapped() ); +} + +void TpetraMatrix::leftRightScale( TpetraVector const & vecLeft, + TpetraVector const & vecRight ) +{ + leftScale( vecLeft ); + rightScale( vecRight ); +} + +void TpetraMatrix::transpose( TpetraMatrix & dst ) const +{ + GEOSX_LAI_ASSERT( ready() ); + + using TpetraTransposer = Tpetra::RowMatrixTransposer< Tpetra_CrsMatrix::scalar_type, + Tpetra_CrsMatrix::local_ordinal_type, + Tpetra_CrsMatrix::global_ordinal_type >; + + TpetraTransposer transposer( Teuchos::RCP< Tpetra_CrsMatrix >( m_matrix.get(), false ) ); + Teuchos::RCP< Tpetra_CrsMatrix > trans = transposer.createTranspose(); + + dst.reset(); + dst.m_matrix.reset( trans.release().get() ); + dst.m_src_map = std::make_unique< Tpetra_Map >( *dst.m_matrix->getDomainMap() ); + dst.m_dst_map = std::make_unique< Tpetra_Map >( *dst.m_matrix->getRangeMap() ); + dst.m_assembled = true; +} + +real64 TpetraMatrix::clearRow( globalIndex const globalRow, + bool const keepDiag, + real64 const diagValue ) +{ + // TODO: Deprecate/remove this method? + + GEOSX_LAI_ASSERT( modifiable() ); + GEOSX_LAI_ASSERT_GE( globalRow, ilower() ); + GEOSX_LAI_ASSERT_GT( iupper(), globalRow ); + + Teuchos::ArrayView< Tpetra_CrsMatrix::local_ordinal_type const > indicesView; + Teuchos::ArrayView< Tpetra_CrsMatrix::scalar_type const > valuesView; + m_matrix->getLocalRowView( getLocalRowID( globalRow ), indicesView, valuesView ); + + array1d< globalIndex > colIndices( indicesView.size() ); + array1d< real64 > values( valuesView.size() ); + + bool const square = numGlobalRows() == numGlobalCols(); + + real64 oldDiag = 0.0; + for( localIndex j = 0; j < indicesView.size(); ++j ) + { + colIndices[j] = m_matrix->getColMap()->getGlobalElement( indicesView[j] ); + if( square && colIndices[j] == globalRow ) + { + oldDiag = valuesView[j]; + values[j] = keepDiag ? oldDiag : diagValue; + } + } + + set( globalRow, colIndices, values ); + return oldDiag; +} + +void TpetraMatrix::addEntries( TpetraMatrix const & src, real64 const scale ) +{ + GEOSX_LAI_ASSERT( ready() ); + GEOSX_LAI_ASSERT( src.ready() ); + GEOSX_LAI_ASSERT( numGlobalRows() == src.numGlobalRows() ); + GEOSX_LAI_ASSERT( numGlobalCols() == src.numGlobalCols() ); + + open(); + Tpetra::MatrixMatrix::Add( src.unwrapped(), false, scale, *m_matrix, 1.0 ); + close(); +} + +void TpetraMatrix::addDiagonal( TpetraVector const & src ) +{ + GEOSX_LAI_ASSERT( ready() ); + GEOSX_LAI_ASSERT( src.ready() ); + GEOSX_LAI_ASSERT( numGlobalRows() == numGlobalCols() ); + GEOSX_LAI_ASSERT( numLocalRows() == src.localSize() ); + + open(); + + // TODO: this is only correct if values are on host + real64 const * const localValues = src.extractLocalVector(); + for( localIndex i = 0; i < numLocalRows(); ++i ) + { + globalIndex const globalRow = getGlobalRowID( i ); + add( globalRow, globalRow, localValues[i] ); + } + + close(); +} + +localIndex TpetraMatrix::maxRowLength() const +{ + GEOSX_LAI_ASSERT( assembled() ); + return m_matrix->getGlobalMaxNumRowEntries(); +} + +localIndex TpetraMatrix::localRowLength( localIndex const localRowIndex ) const +{ + GEOSX_LAI_ASSERT( assembled() ); + return m_matrix->getNumEntriesInLocalRow( localRowIndex ); +} + +localIndex TpetraMatrix::globalRowLength( globalIndex const globalRowIndex ) const +{ + GEOSX_LAI_ASSERT( assembled() ); + return m_matrix->getNumEntriesInGlobalRow( globalRowIndex ); +} + +void TpetraMatrix::getRowCopy( globalIndex const globalRow, + arraySlice1d< globalIndex > const & colIndices, + arraySlice1d< real64 > const & values ) const +{ + // TODO: Deprecate/remove this method? + + GEOSX_LAI_ASSERT( ready() ); + GEOSX_LAI_ASSERT_GE( globalRow, ilower() ); + GEOSX_LAI_ASSERT_GT( iupper(), globalRow ); + + Teuchos::ArrayView< Tpetra_CrsMatrix::local_ordinal_type const > indicesView; + Teuchos::ArrayView< Tpetra_CrsMatrix::scalar_type const > valuesView; + m_matrix->getLocalRowView( getLocalRowID( globalRow ), indicesView, valuesView ); + + GEOSX_LAI_ASSERT_GE( colIndices.size(), indicesView.size() ); + GEOSX_LAI_ASSERT_GE( values.size(), indicesView.size() ); + + localIndex const numEntries = LvArray::integerConversion< localIndex >( indicesView.size() ); + for( localIndex i = 0; i < numEntries; ++i ) + { + colIndices[i] = LvArray::integerConversion< globalIndex >( m_src_map->getGlobalElement( indicesView[i] ) ); + values[i] = valuesView[i]; + } +} + +real64 TpetraMatrix::getDiagValue( globalIndex const globalRow ) const +{ + // TODO: Deprecate/remove this method? + + GEOSX_LAI_ASSERT( assembled() ); + GEOSX_LAI_ASSERT_GE( globalRow, ilower() ); + GEOSX_LAI_ASSERT_GT( iupper(), globalRow ); + + Teuchos::ArrayView< Tpetra_CrsMatrix::local_ordinal_type const > indicesView; + Teuchos::ArrayView< Tpetra_CrsMatrix::scalar_type const > valuesView; + m_matrix->getLocalRowView( getLocalRowID( globalRow ), indicesView, valuesView ); + + localIndex const numEntries = LvArray::integerConversion< localIndex >( indicesView.size() ); + for( localIndex j = 0; j < numEntries; ++j ) + { + if( m_src_map->getGlobalElement( indicesView[j] ) == globalRow ) + { + return valuesView[j]; + } + } + + return 0.0; +} + +void TpetraMatrix::extractDiagonal( TpetraVector & dst ) const +{ + GEOSX_LAI_ASSERT( ready() ); + GEOSX_LAI_ASSERT( dst.ready() ); + GEOSX_LAI_ASSERT_EQ( dst.localSize(), numLocalRows() ); + + m_matrix->getLocalDiagCopy( dst.unwrapped() ); +} + +TpetraMatrix::Tpetra_CrsMatrix const & TpetraMatrix::unwrapped() const +{ + GEOSX_LAI_ASSERT( created() ); + return *m_matrix; +} + +TpetraMatrix::Tpetra_CrsMatrix & TpetraMatrix::unwrapped() +{ + GEOSX_LAI_ASSERT( created() ); + return *m_matrix; +} + +globalIndex TpetraMatrix::numGlobalRows() const +{ + GEOSX_LAI_ASSERT( created() ); + return LvArray::integerConversion< globalIndex >( m_dst_map->getGlobalNumElements() ); +} + +globalIndex TpetraMatrix::numGlobalCols() const +{ + GEOSX_LAI_ASSERT( created() ); + return LvArray::integerConversion< globalIndex >( m_src_map->getGlobalNumElements() ); +} + +globalIndex TpetraMatrix::ilower() const +{ + GEOSX_LAI_ASSERT( created() ); + return m_dst_map->getMinGlobalIndex(); +} + +globalIndex TpetraMatrix::iupper() const +{ + GEOSX_LAI_ASSERT( created() ); + return m_dst_map->getMaxGlobalIndex() + 1; +} + +globalIndex TpetraMatrix::jlower() const +{ + GEOSX_LAI_ASSERT( created() ); + return m_src_map->getMinGlobalIndex(); +} + +globalIndex TpetraMatrix::jupper() const +{ + GEOSX_LAI_ASSERT( created() ); + return m_src_map->getMaxGlobalIndex() + 1; +} + +localIndex TpetraMatrix::numLocalNonzeros() const +{ + GEOSX_LAI_ASSERT( assembled() ); + return LvArray::integerConversion< localIndex >( m_matrix->getNodeNumEntries() ); +} + +globalIndex TpetraMatrix::numGlobalNonzeros() const +{ + GEOSX_LAI_ASSERT( assembled() ); + return LvArray::integerConversion< globalIndex >( m_matrix->getGlobalNumEntries() ); +} + +real64 TpetraMatrix::normInf() const +{ + GEOSX_LAI_ASSERT( ready() ); + + using LocalMatrix = Tpetra_CrsMatrix::local_matrix_type; + LocalMatrix const & localMatrix = m_matrix->getLocalMatrix(); + + real64 localNorm = 0.0; + Kokkos::parallel_reduce( Kokkos::RangePolicy< Tpetra_CrsMatrix::execution_space >( 0, localMatrix.numRows() ), + KOKKOS_LAMBDA ( localIndex const i, real64 & maxVal ) + { + KokkosSparse::SparseRowViewConst< LocalMatrix > const rowView = localMatrix.rowConst( i ); + real64 rowSumAbs = 0.0; + for( int k = 0; k < rowView.length; ++k ) + { + rowSumAbs += fabs( rowView.value( k ) ); + } + maxVal = fmax( maxVal, rowSumAbs ); + }, Kokkos::Max< real64, Kokkos::HostSpace >( localNorm ) ); + + return MpiWrapper::Max( localNorm, getComm() ); +} + +real64 TpetraMatrix::norm1() const +{ + GEOSX_LAI_ASSERT( ready() ); + + TpetraMatrix matTrans; + transpose( matTrans ); + return matTrans.normInf(); +} + +real64 TpetraMatrix::normFrobenius() const +{ + GEOSX_LAI_ASSERT( ready() ); + return m_matrix->getFrobeniusNorm(); +} + +localIndex TpetraMatrix::getLocalRowID( globalIndex const globalRow ) const +{ + GEOSX_LAI_ASSERT( created() ); + return m_dst_map->getLocalElement( globalRow ); +} + +globalIndex TpetraMatrix::getGlobalRowID( localIndex const localRow ) const +{ + GEOSX_LAI_ASSERT( created() ); + GEOSX_LAI_ASSERT_GE( localRow, 0 ); + GEOSX_LAI_ASSERT_GT( numLocalRows(), localRow ); + return m_dst_map->getGlobalElement( localRow ); +} + +localIndex TpetraMatrix::numLocalCols() const +{ + GEOSX_LAI_ASSERT( created() ); + return LvArray::integerConversion< localIndex >( m_src_map->getNodeNumElements() ); +} + +localIndex TpetraMatrix::numLocalRows() const +{ + GEOSX_LAI_ASSERT( created() ); + return LvArray::integerConversion< localIndex >( m_dst_map->getNodeNumElements() ); +} + +MPI_Comm TpetraMatrix::getComm() const +{ + GEOSX_LAI_ASSERT( created() ); +#ifdef GEOSX_USE_MPI + return *dynamic_cast< Tpetra_Comm const & >( *m_matrix->getRowMap()->getComm() ).getRawMpiComm(); +#else + return MPI_COMM_GEOSX; +#endif +} + +void TpetraMatrix::print( std::ostream & os ) const +{ + GEOSX_LAI_ASSERT( ready() ); + Teuchos::RCP< Teuchos::FancyOStream > stream = Teuchos::getFancyOStream( Teuchos::rcp( &os, false ) ); + m_matrix->describe( *stream, Teuchos::VERB_EXTREME ); +} + +void TpetraMatrix::write( string const & filename, + LAIOutputFormat const format ) const +{ + GEOSX_LAI_ASSERT( ready() ); + switch( format ) + { + case LAIOutputFormat::NATIVE_ASCII: + { + std::ofstream ofs( filename ); + print( ofs ); + break; + } + case LAIOutputFormat::MATRIX_MARKET: + { + using Writer = Tpetra::MatrixMarket::Writer< Tpetra_CrsMatrix >; + Writer::writeSparseFile( filename, *m_matrix ); + break; + } + default: + { + GEOSX_ERROR( "Unsupported matrix output format" ); + } + } +} + +void TpetraMatrix::multiply( bool const transA, + TpetraMatrix const & B, + bool const transB, + TpetraMatrix & C ) const +{ + GEOSX_LAI_ASSERT( ready() ); + GEOSX_LAI_ASSERT( B.ready() ); + + C.createWithLocalSize( transA ? numLocalCols() : numLocalRows(), + transB ? B.numLocalRows() : B.numLocalCols(), + 1, // TODO: estimate entries per row? + getComm() ); + + Tpetra::MatrixMatrix::Multiply( unwrapped(), + transA, + B.unwrapped(), + transB, + C.unwrapped(), + true ); + C.m_assembled = true; +} + +} // end geosx namespace diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/TpetraMatrix.hpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TpetraMatrix.hpp new file mode 100644 index 00000000000..4c5597a1c76 --- /dev/null +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TpetraMatrix.hpp @@ -0,0 +1,328 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2019 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2019 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2019 Total, S.A + * Copyright (c) 2019- GEOSX Contributors + * All right reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file TpetraMatrix.hpp + */ + +#ifndef GEOSX_LINEARALGEBRA_INTERFACES_TPETRAMATRIX_HPP +#define GEOSX_LINEARALGEBRA_INTERFACES_TPETRAMATRIX_HPP + +#include "linearAlgebra/interfaces/LinearOperator.hpp" +#include "linearAlgebra/interfaces/MatrixBase.hpp" +#include "linearAlgebra/interfaces/trilinos/TpetraVector.hpp" + +#include +#include + +namespace geosx +{ + +/** + * @brief Wrapper class for Trilinos/Tpetra's CrsMatrix class. + */ +class TpetraMatrix : public virtual LinearOperator< TpetraVector >, + private MatrixBase< TpetraMatrix, TpetraVector > +{ +public: + + /// Compatible vector type + using Vector = TpetraVector; + + /** + * @name Constructor/Destructor methods + */ + ///@{ + + /** + * @brief Empty matrix constructor. + * + * Create an empty (distributed) matrix. + */ + TpetraMatrix(); + + /** + * @brief Copy constructor. + * @param[in] src the matrix to be copied + */ + TpetraMatrix( TpetraMatrix const & src ); + + /** + * @brief Virtual destructor. + */ + virtual ~TpetraMatrix() override; + + ///@} + + /** + * @name MatrixBase interface + */ + ///@{ + + using MatrixBase::createWithLocalSize; + using MatrixBase::createWithGlobalSize; + using MatrixBase::create; + using MatrixBase::closed; + using MatrixBase::assembled; + using MatrixBase::insertable; + using MatrixBase::modifiable; + using MatrixBase::ready; + using MatrixBase::residual; + + virtual void createWithLocalSize( localIndex const localRows, + localIndex const localCols, + localIndex const maxEntriesPerRow, + MPI_Comm const & comm ) override; + + virtual void createWithGlobalSize( globalIndex const globalRows, + globalIndex const globalCols, + localIndex const maxEntriesPerRow, + MPI_Comm const & comm ) override; + + virtual void open() override; + + virtual void close() override; + + virtual bool created() const override; + + virtual void reset() override; + + virtual void set( real64 const value ) override; + + virtual void zero() override; + + virtual void add( globalIndex const rowIndex, + globalIndex const colIndex, + real64 const value ) override; + + virtual void set( globalIndex const rowIndex, + globalIndex const colIndex, + real64 const value ) override; + + virtual void insert( globalIndex const rowIndex, + globalIndex const colIndex, + real64 const value ) override; + + virtual void add( globalIndex const rowIndex, + globalIndex const * colIndices, + real64 const * values, + localIndex const size ) override; + + virtual void set( globalIndex const rowIndex, + globalIndex const * colIndices, + real64 const * values, + localIndex const size ) override; + + virtual void insert( globalIndex const rowIndex, + globalIndex const * colIndices, + real64 const * values, + localIndex const size ) override; + + virtual void add( globalIndex const rowIndex, + arraySlice1d< globalIndex const > const & colIndices, + arraySlice1d< real64 const > const & values ) override; + + virtual void set( globalIndex const rowIndex, + arraySlice1d< globalIndex const > const & colIndices, + arraySlice1d< real64 const > const & values ) override; + + virtual void insert( globalIndex const rowIndex, + arraySlice1d< globalIndex const > const & colIndices, + arraySlice1d< real64 const > const & values ) override; + + virtual void add( arraySlice1d< globalIndex const > const & rowIndices, + arraySlice1d< globalIndex const > const & colIndices, + arraySlice2d< real64 const, MatrixLayout::ROW_MAJOR > const & values ) override; + + virtual void set( arraySlice1d< globalIndex const > const & rowIndices, + arraySlice1d< globalIndex const > const & colIndices, + arraySlice2d< real64 const, MatrixLayout::ROW_MAJOR > const & values ) override; + + virtual void insert( arraySlice1d< globalIndex const > const & rowIndices, + arraySlice1d< globalIndex const > const & colIndices, + arraySlice2d< real64 const, MatrixLayout::ROW_MAJOR > const & values ) override; + + virtual void add( globalIndex const * rowIndices, + globalIndex const * colIndices, + real64 const * values, + localIndex const numRows, + localIndex const numCols ) override; + + virtual void set( globalIndex const * rowIndices, + globalIndex const * colIndices, + real64 const * values, + localIndex const numRows, + localIndex const numCols ) override; + + virtual void insert( globalIndex const * rowIndices, + globalIndex const * colIndices, + real64 const * values, + localIndex const numRows, + localIndex const numCols ) override; + + virtual void apply( TpetraVector const & src, + TpetraVector & dst ) const override; + + virtual void applyTranspose( TpetraVector const & src, + TpetraVector & dst ) const override; + + virtual void multiply( TpetraMatrix const & src, + TpetraMatrix & dst ) const override; + + virtual void leftMultiplyTranspose( TpetraMatrix const & src, + TpetraMatrix & dst ) const override; + + virtual void rightMultiplyTranspose( TpetraMatrix const & src, + TpetraMatrix & dst ) const override; + + virtual void multiplyRAP( TpetraMatrix const & R, + TpetraMatrix const & P, + TpetraMatrix & dst ) const override; + + virtual void multiplyPtAP( TpetraMatrix const & P, + TpetraMatrix & dst ) const override; + + virtual void gemv( real64 const alpha, + TpetraVector const & x, + real64 const beta, + TpetraVector & y, + bool useTranspose = false ) const override; + + virtual void scale( real64 const scalingFactor ) override; + + virtual void leftScale( TpetraVector const & vec ) override; + + virtual void rightScale( TpetraVector const & vec ) override; + + virtual void leftRightScale( TpetraVector const & vecLeft, + TpetraVector const & vecRight ) override; + + virtual void transpose( TpetraMatrix & dst ) const override; + + virtual real64 clearRow( globalIndex const row, + bool const keepDiag = false, + real64 const diagValue = 0.0 ) override; + + virtual void addEntries( TpetraMatrix const & src, real64 const scale = 1.0 ) override; + + virtual void addDiagonal( TpetraVector const & src ) override; + + virtual localIndex maxRowLength() const override; + + virtual localIndex localRowLength( localIndex localRowIndex ) const override; + + virtual localIndex globalRowLength( globalIndex globalRowIndex ) const override; + + virtual void getRowCopy( globalIndex globalRow, + arraySlice1d< globalIndex > const & colIndices, + arraySlice1d< real64 > const & values ) const override; + + virtual real64 getDiagValue( globalIndex globalRow ) const override; + + virtual void extractDiagonal( TpetraVector & dst ) const override; + + virtual globalIndex numGlobalRows() const override; + + virtual globalIndex numGlobalCols() const override; + + virtual localIndex numLocalRows() const override; + + virtual localIndex numLocalCols() const override; + + virtual globalIndex ilower() const override; + + virtual globalIndex iupper() const override; + + virtual globalIndex jlower() const override; + + virtual globalIndex jupper() const override; + + virtual localIndex numLocalNonzeros() const override; + + virtual globalIndex numGlobalNonzeros() const override; + + virtual real64 normInf() const override; + + virtual real64 norm1() const override; + + virtual real64 normFrobenius() const override; + + virtual localIndex getLocalRowID( globalIndex const globalRow ) const override; + + virtual globalIndex getGlobalRowID( localIndex const localRow ) const override; + + virtual MPI_Comm getComm() const override; + + virtual void print( std::ostream & os = std::cout ) const override; + + virtual void write( string const & filename, + LAIOutputFormat const format = LAIOutputFormat::MATRIX_MARKET ) const override; + + ///@} + + /** + * @brief Alias for Tpetra map template instantiation used by this class. + */ + using Tpetra_Map = Tpetra::Map< int, globalIndex >; + + /** + * @brief Alias for specific Tpetra matrix template instantiation wrapped by this class. + * + * @note This uses Tpetra's default execution/memory space. When built with CUDA support, + * this will be equal to Kokkos::Cuda, so we won't be able to create a host-only vector. + * If we want both in the same executable, we'll have to make adjustments to our LAI approach. + */ + using Tpetra_CrsMatrix = Tpetra::CrsMatrix< real64, int, globalIndex >; + + /** + * @brief Returns a const pointer to the underlying matrix. + * @return const pointer to the underlying matrix + */ + Tpetra_CrsMatrix const & unwrapped() const; + + /** + * @brief Returns a non-const pointer to the underlying matrix. + * @return non-const pointer to the underlying matrix + */ + Tpetra_CrsMatrix & unwrapped(); + +private: + + /** + * @brief Perform a matrix matrix product with Parallel Matrix + */ + void multiply( bool const transA, + TpetraMatrix const & B, + bool const transB, + TpetraMatrix & C ) const; + + /** + * @brief Create the matrix by copying data from an Epetra_CrsMatrix + * @param src the source matrix + */ + void create( Tpetra_CrsMatrix const & src ); + + /// Pointer to the underlying Epetra_CrsMatrix. + std::unique_ptr< Tpetra_CrsMatrix > m_matrix; + + /// Map representing the parallel partitioning of a source vector (x in y=Ax) + std::unique_ptr< Tpetra_Map > m_src_map; + + /// Map representing the parallel partitioning of a destination vector (y in y=Ax) + std::unique_ptr< Tpetra_Map > m_dst_map; +}; + +} // namespace geosx + +#endif //GEOSX_LINEARALGEBRA_INTERFACES_TPETRAMATRIX_HPP diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/TpetraUtils.hpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TpetraUtils.hpp new file mode 100644 index 00000000000..832697fbf1d --- /dev/null +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TpetraUtils.hpp @@ -0,0 +1,35 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2019 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2019 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2019 Total, S.A + * Copyright (c) 2019- GEOSX Contributors + * All right reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file TpetraUtils.hpp + */ +#ifndef GEOSX_TPETRAUTILS_HPP +#define GEOSX_TPETRAUTILS_HPP + +#include "common/GeosxConfig.hpp" + +#ifdef GEOSX_USE_MPI +#include + +/// Alias for specific Tpetra::Comm implementation used. +using Tpetra_Comm = Teuchos::MpiComm< int >; +#else +#include + +/// Alias for specific Tpetra::Comm implementation used. +using Tpetra_Comm = Teuchos::SerialComm< int >; +#endif + +#endif //GEOSX_TPETRAUTILS_HPP diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/TpetraVector.cpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TpetraVector.cpp new file mode 100644 index 00000000000..7585892e072 --- /dev/null +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TpetraVector.cpp @@ -0,0 +1,483 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2019 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2019 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2019 Total, S.A + * Copyright (c) 2019- GEOSX Contributors + * All right reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file TpetraVector.cpp + */ + +#include "TpetraVector.hpp" + +#include "codingUtilities/Utilities.hpp" +#include "linearAlgebra/interfaces/trilinos/TpetraUtils.hpp" + +#include +#include +#include + +namespace geosx +{ + +TpetraVector::TpetraVector() + : VectorBase(), + m_vector{} +{} + +TpetraVector::TpetraVector( TpetraVector const & src ) + : TpetraVector() +{ + *this = src; +} + +TpetraVector::TpetraVector( TpetraVector && src ) noexcept + : TpetraVector() +{ + *this = std::move( src ); +} + +TpetraVector & TpetraVector::operator=( TpetraVector const & src ) +{ + GEOSX_LAI_ASSERT( &src != this ); + GEOSX_LAI_ASSERT( src.ready() ); + if( m_vector ) + { + Tpetra::deep_copy( *m_vector, *src.m_vector ); + } + else + { + m_vector = std::make_unique< Tpetra_Vector >( *src.m_vector, Teuchos::Copy ); + } + return *this; +} + +TpetraVector & TpetraVector::operator=( TpetraVector && src ) noexcept +{ + GEOSX_LAI_ASSERT( &src != this ); + GEOSX_LAI_ASSERT( src.ready() ); + m_vector = std::move( src.m_vector ); + return *this; +} + +TpetraVector::~TpetraVector() = default; + +bool TpetraVector::created() const +{ + return bool(m_vector); +} + +void TpetraVector::createWithLocalSize( localIndex const localSize, + MPI_Comm const & MPI_PARAM( comm ) ) +{ + GEOSX_LAI_ASSERT( closed() ); + GEOSX_LAI_ASSERT_GE( localSize, 0 ); + + Teuchos::RCP< Tpetra_Comm > tpetraComm( new Tpetra_Comm( MPI_PARAM( comm ) ) ); + Teuchos::RCP< Tpetra_Map > map( new Tpetra_Map( Teuchos::OrdinalTraits< Tpetra::global_size_t >::invalid(), + LvArray::integerConversion< size_t >( localSize ), + 0, + tpetraComm ) ); + m_vector = std::make_unique< Tpetra_Vector >( map ); +} + +void TpetraVector::createWithGlobalSize( globalIndex const globalSize, + MPI_Comm const & MPI_PARAM( comm ) ) +{ + GEOSX_LAI_ASSERT( closed() ); + GEOSX_LAI_ASSERT_GE( globalSize, 0 ); + + Teuchos::RCP< Tpetra_Comm > tpetraComm( new Tpetra_Comm( MPI_PARAM( comm ) ) ); + Teuchos::RCP< Tpetra_Map > map( new Tpetra_Map( LvArray::integerConversion< Tpetra::global_size_t >( globalSize ), + 0, + tpetraComm ) ); + m_vector = std::make_unique< Tpetra_Vector >( map ); +} + +void TpetraVector::createWithLocalValues( arrayView1d< real64 > const & localValues, + MPI_Comm const & comm ) +{ + Teuchos::RCP< Tpetra_Comm > tpetraComm( new Tpetra_Comm( MPI_PARAM( comm ) ) ); + Teuchos::RCP< Tpetra_Map > map( new Tpetra_Map( Teuchos::OrdinalTraits< Tpetra::global_size_t >::invalid(), + LvArray::integerConversion< size_t >( localValues.size() ), + 0, + tpetraComm ) ); + + // Define what "device" is +#ifdef GEOSX_USE_CUDA + using KokkosSpaceDevice = Kokkos::CudaSpace; // Not Kokkos::CudaUVMSpace, since we wrap a normal device pointer + localValues.move( LvArray::MemorySpace::GPU, false ); // Make sure data is current on device +#else + using KokkosSpaceDevice = Kokkos::HostSpace; +#endif // GEOSX_USE_CUDA + + // Wrap a non-owning Kokkos view around input data + using InputView = Kokkos::View< real64 * *, Kokkos::LayoutLeft, KokkosSpaceDevice, Kokkos::MemoryUnmanaged >; + InputView const localValuesView( localValues.data(), localValues.size(), 1 ); + + // Sadly, Tpetra currently **requires** the use of UM as default memory space for all its views. + // The disabled code below would be valid if we were allowed to turn UM off and use regular device allocation. + // Unfortunately, we are not, and zero-copy construction from a device pointer is impossible to implement. + // TODO: reevaluate this if Tpetra relaxes the UM requirement. +#if 0 + + // Specify a non-owning DualView + using KokkosDualView = Kokkos::DualView< real64 * *, // this type required by Tpetra::MultiVector, even for a 1D vector + Kokkos::LayoutLeft, // doesn't matter, since the first extent will be 1 + Tpetra_Vector::device_type::execution_space, + Kokkos::MemoryUnmanaged >; // tell the view we own the memory + + // DualView constructor requires that host and device views are in sync. + localValues.move( LvArray::MemorySpace::CPU, false ); + + // Directly use device view of input data (shallow copy) + KokkosDualView::t_dev const values_d = localValuesView; + + // Also create a non-owning view pointing to host data + KokkosDualView::t_host const values_h( localValues.data( LvArray::MemorySpace::CPU ), localValues.size(), 1 ); + +#else + + // Use an owning view, since we'll be making a copy + using KokkosDualView = Tpetra_Vector::dual_view_type; + + // Make a deep copy of input values on device + KokkosDualView::t_dev const values_d( "", localValuesView.size(), 1 ); + Kokkos::deep_copy( values_d, localValuesView ); + + // Also make a host copy as required by DualView + KokkosDualView::t_host const values_h = Kokkos::create_mirror_view( values_d ); + +#endif // 0 + + // Make a DualView containing host and device views (either shallow or deep) and feed to vector + KokkosDualView values( values_d, values_h ); + m_vector = std::make_unique< Tpetra_Vector >( map, values ); +} + +void TpetraVector::open() +{ + GEOSX_LAI_ASSERT( ready() ); + m_closed = false; + m_vector->sync_host(); // since all element-wise modifications work on host +} + +void TpetraVector::close() +{ + GEOSX_LAI_ASSERT( !closed() ); + m_closed = true; + m_vector->sync_device(); // assuming we've modified host data while vector was open + + // Nothing to do w.r.t. parallel sync since Tpetra::Vector does not support remote assembly. + // Tpetra::FEMultiVector does, but as of May 2020 it is still lacking proper constructors, documentation, etc. +} + +void TpetraVector::reset() +{ + VectorBase::reset(); + m_vector.reset(); +} + +void TpetraVector::set( globalIndex const globalRowIndex, + real64 const value ) +{ + GEOSX_LAI_ASSERT( !closed() ); + m_vector->replaceGlobalValue( globalRowIndex, value ); + m_vector->modify_host(); +} + +void TpetraVector::add( globalIndex const globalRowIndex, + real64 const value ) +{ + GEOSX_LAI_ASSERT( !closed() ); + m_vector->sumIntoGlobalValue( globalRowIndex, value ); + m_vector->modify_host(); +} + +void TpetraVector::add( globalIndex const * globalRowIndices, + real64 const * values, + localIndex const size ) +{ + GEOSX_LAI_ASSERT( !closed() ); + for( localIndex i = 0; i < size; ++i ) + { + m_vector->sumIntoGlobalValue( globalRowIndices[i], values[i] ); + } + m_vector->modify_host(); +} + +void TpetraVector::set( globalIndex const * globalRowIndices, + real64 const * values, + localIndex const size ) +{ + GEOSX_LAI_ASSERT( !closed() ); + for( localIndex i = 0; i < size; ++i ) + { + m_vector->replaceGlobalValue( globalRowIndices[i], values[i] ); + } + m_vector->modify_host(); +} + +void TpetraVector::set( arraySlice1d< globalIndex const > const & globalRowIndices, + arraySlice1d< real64 const > const & values ) +{ + set( globalRowIndices.dataIfContiguous(), values.dataIfContiguous(), globalRowIndices.size() ); +} + +void TpetraVector::add( arraySlice1d< globalIndex const > const & globalRowIndices, + arraySlice1d< real64 const > const & values ) +{ + add( globalRowIndices.dataIfContiguous(), values.dataIfContiguous(), globalRowIndices.size() ); +} + +void TpetraVector::set( real64 const value ) +{ + GEOSX_LAI_ASSERT( ready() ); + m_vector->modify_device(); // ensure putScalar runs on device + m_vector->putScalar( value ); + m_vector->sync_host(); +} + +void TpetraVector::zero() +{ + set( 0.0 ); +} + +void TpetraVector::rand( unsigned int const GEOSX_UNUSED_PARAM( seed ) ) +{ + GEOSX_LAI_ASSERT( ready() ); + // There's no way currently (May 2020) to set the random seed. + m_vector->randomize( -1.0, 1.0 ); +} + +void TpetraVector::scale( real64 const scalingFactor ) +{ + GEOSX_LAI_ASSERT( ready() ); + m_vector->scale( scalingFactor ); +} + +void TpetraVector::reciprocal() +{ + GEOSX_LAI_ASSERT( ready() ); + m_vector->reciprocal( *m_vector ); +} + +real64 TpetraVector::dot( TpetraVector const & vec ) const +{ + GEOSX_LAI_ASSERT( ready() ); + GEOSX_LAI_ASSERT( vec.ready() ); + GEOSX_LAI_ASSERT_EQ( globalSize(), vec.globalSize() ); + + return m_vector->dot( *vec.m_vector ); +} + +void TpetraVector::copy( TpetraVector const & x ) +{ + GEOSX_LAI_ASSERT( ready() ); + GEOSX_LAI_ASSERT( x.ready() ); + GEOSX_LAI_ASSERT_EQ( globalSize(), x.globalSize() ); + + Tpetra::deep_copy( *m_vector, *x.m_vector ); +} + +void TpetraVector::axpy( real64 alpha, + TpetraVector const & x ) +{ + GEOSX_LAI_ASSERT( ready() ); + GEOSX_LAI_ASSERT( x.ready() ); + GEOSX_LAI_ASSERT_EQ( globalSize(), x.globalSize() ); + + m_vector->update( alpha, x.unwrapped(), 1. ); +} + +void TpetraVector::axpby( real64 alpha, + TpetraVector const & x, + real64 beta ) +{ + GEOSX_LAI_ASSERT( ready() ); + GEOSX_LAI_ASSERT( x.ready() ); + GEOSX_LAI_ASSERT_EQ( globalSize(), x.globalSize() ); + + m_vector->update( alpha, x.unwrapped(), beta ); +} + +real64 TpetraVector::norm1() const +{ + GEOSX_LAI_ASSERT( ready() ); + return m_vector->norm1(); +} + +real64 TpetraVector::norm2() const +{ + GEOSX_LAI_ASSERT( ready() ); + return m_vector->norm2(); +} + +real64 TpetraVector::normInf() const +{ + GEOSX_LAI_ASSERT( ready() ); + return m_vector->normInf(); +} + +globalIndex TpetraVector::globalSize() const +{ + GEOSX_LAI_ASSERT( created() ); + return m_vector->getGlobalLength(); +} + +localIndex TpetraVector::localSize() const +{ + GEOSX_LAI_ASSERT( created() ); + return m_vector->getLocalLength(); +} + +globalIndex TpetraVector::ilower() const +{ + GEOSX_LAI_ASSERT( created() ); + return m_vector->getMap()->getMinGlobalIndex(); +} + +globalIndex TpetraVector::iupper() const +{ + GEOSX_LAI_ASSERT( created() ); + return m_vector->getMap()->getMaxGlobalIndex() + 1; +} + +real64 TpetraVector::get( globalIndex const globalRow ) const +{ + GEOSX_LAI_ASSERT( ready() ); + GEOSX_LAI_ASSERT_GE( globalRow, ilower() ); + GEOSX_LAI_ASSERT_GT( iupper(), globalRow ); + + // TODO: deprecate? + + m_vector->sync_host(); + return m_vector->getLocalViewHost().data()[ getLocalRowID( globalRow ) ]; +} + +void TpetraVector::get( arraySlice1d< globalIndex const > const & globalIndices, + arraySlice1d< real64 > const & values ) const +{ + GEOSX_LAI_ASSERT( ready() ); + GEOSX_LAI_ASSERT_GE( values.size(), globalIndices.size() ); + GEOSX_LAI_ASSERT_GE( *std::min_element( globalIndices.dataIfContiguous(), globalIndices.dataIfContiguous() + globalIndices.size() ), ilower() ); + GEOSX_LAI_ASSERT_GT( iupper(), *std::max_element( globalIndices.dataIfContiguous(), + globalIndices.dataIfContiguous() + globalIndices.size() ) ); + + // TODO: deprecate? + + m_vector->sync_host(); + real64 const * const localVector = m_vector->getLocalViewHost().data(); + for( localIndex i = 0; i < globalIndices.size(); ++i ) + { + values[i] = localVector[ getLocalRowID( globalIndices[i] ) ]; + } +} + +localIndex TpetraVector::getLocalRowID( globalIndex const globalRow ) const +{ + GEOSX_LAI_ASSERT( created() ); + return m_vector->getMap()->getLocalElement( globalRow ); +} + +globalIndex TpetraVector::getGlobalRowID( localIndex const localRow ) const +{ + GEOSX_LAI_ASSERT( created() ); + return m_vector->getMap()->getGlobalElement( localRow ); +} + +real64 const * TpetraVector::extractLocalVector() const +{ + GEOSX_LAI_ASSERT( ready() ); + + // TODO: deprecate? + + // TODO: This will always return pointer in vector's target memory space (i.e. device with CUDA). + // If we want to get data in a different space, we'll need separate methods or a parameter. + return m_vector->getLocalViewDevice().data(); +} + +real64 * TpetraVector::extractLocalVector() +{ + GEOSX_LAI_ASSERT( ready() ); + + // TODO: deprecate? + + // TODO: This will always return pointer in vector's target memory space (i.e. device with CUDA). + // If we want to get data in a different space, we'll need separate methods or a parameter. + + // Mark data dirty since we are returning a pointer to non-const to the user + m_vector->modify_device(); + return m_vector->getLocalViewDevice().data(); +} + +void TpetraVector::extract( arrayView1d< real64 > const & localVector ) const +{ + GEOSX_LAI_ASSERT_EQ( localSize(), localVector.size() ); + m_vector->sync_device(); + real64 const * const data = m_vector->getLocalViewDevice().data(); + forAll< parallelDevicePolicy<> >( localSize(), [=] GEOSX_DEVICE ( localIndex const k ) + { + localVector[k] = data[k]; + } ); +} + +MPI_Comm TpetraVector::getComm() const +{ + GEOSX_LAI_ASSERT( created() ); +#ifdef GEOSX_USE_MPI + return *dynamic_cast< Tpetra_Comm const & >( *m_vector->getMap()->getComm() ).getRawMpiComm(); +#else + return MPI_COMM_GEOSX; +#endif +} + +void TpetraVector::print( std::ostream & os ) const +{ + GEOSX_LAI_ASSERT( ready() ); + Teuchos::RCP< Teuchos::FancyOStream > stream = Teuchos::getFancyOStream( Teuchos::rcp( &os, false ) ); + m_vector->describe( *stream, Teuchos::VERB_EXTREME ); +} + +void TpetraVector::write( string const & filename, + LAIOutputFormat format ) const +{ + GEOSX_LAI_ASSERT( ready() ); + switch( format ) + { + case LAIOutputFormat::MATRIX_MARKET: + { + using Tpetra_CrsMatrix = Tpetra::CrsMatrix< Tpetra_Vector::scalar_type, Tpetra_Vector::local_ordinal_type, Tpetra_Vector::global_ordinal_type >; + using Writer = Tpetra::MatrixMarket::Writer< Tpetra_CrsMatrix >; + Writer::writeDenseFile( filename, *m_vector ); + break; + } + default: + { + GEOSX_ERROR( "Unsupported vector output format" ); + } + } +} + +TpetraVector::Tpetra_Vector const & TpetraVector::unwrapped() const +{ + GEOSX_LAI_ASSERT( created() ); + return *m_vector; +} + +TpetraVector::Tpetra_Vector & TpetraVector::unwrapped() +{ + GEOSX_LAI_ASSERT( created() ); + return *m_vector; +} + + +} // namespace geosx diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/TpetraVector.hpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TpetraVector.hpp new file mode 100644 index 00000000000..19839467431 --- /dev/null +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TpetraVector.hpp @@ -0,0 +1,228 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2019 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2019 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2019 Total, S.A + * Copyright (c) 2019- GEOSX Contributors + * All right reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file TpetraVector.hpp + */ + +#ifndef GEOSX_LINEARALGEBRA_INTERFACES_TPETRAVECTOR_HPP_ +#define GEOSX_LINEARALGEBRA_INTERFACES_TPETRAVECTOR_HPP_ + +#include "linearAlgebra/interfaces/VectorBase.hpp" + +#include +#include +#include + +namespace geosx +{ + +/** + * @brief Wrapper class for Trilinos/Tpetra's Vector class. + */ +class TpetraVector final : private VectorBase< TpetraVector > +{ +public: + + /** + * @name Constructor/Destructor Methods + */ + ///@{ + + /** + * @brief Empty vector constructor. + * Create an empty (distributed) vector. + */ + TpetraVector(); + + /** + * @brief Copy constructor. + * @param src vector to be copied + */ + TpetraVector( TpetraVector const & src ); + + /** + * @brief Move constructor + * @param src vector to move from + */ + TpetraVector( TpetraVector && src ) noexcept; + + /** + * @brief Copy assignment. + * @param src vector to be copied + * @return reference to this object + */ + TpetraVector & operator=( TpetraVector const & src ); + + /** + * @brief Move assignment. + * @param src vector to move from + * @return reference to this object + */ + TpetraVector & operator=( TpetraVector && src ) noexcept; + + /** + * @brief Destructor. + */ + ~TpetraVector(); + + ///@} + + /** + * @name VectorBase interface + */ + ///@{ + + using VectorBase::closed; + using VectorBase::ready; + + virtual bool created() const override; + + virtual void createWithLocalSize( localIndex const localSize, + MPI_Comm const & comm ) override; + + virtual void createWithGlobalSize( globalIndex const globalSize, + MPI_Comm const & comm ) override; + + virtual void createWithLocalValues( arrayView1d< real64 > const & localValues, + MPI_Comm const & comm ) override; + + virtual void open() override; + + virtual void close() override; + + virtual void reset() override; + + virtual void set( globalIndex const globalRowIndex, + real64 const value ) override; + + virtual void add( globalIndex const globalRowIndex, + real64 const value ) override; + + virtual void set( globalIndex const * globalRowIndices, + real64 const * values, + localIndex size ) override; + + virtual void add( globalIndex const * globalRowIndices, + real64 const * values, + localIndex const size ) override; + + virtual void set( arraySlice1d< globalIndex const > const & globalRowIndices, + arraySlice1d< real64 const > const & values ) override; + + virtual void add( arraySlice1d< globalIndex const > const & globalRowIndices, + arraySlice1d< real64 const > const & values ) override; + + virtual void set( real64 const value ) override; + + virtual void zero() override; + + virtual void rand( unsigned const seed = 1984 ) override; + + virtual void scale( real64 const scalingFactor ) override; + + virtual void reciprocal() override; + + virtual real64 dot( TpetraVector const & vec ) const override; + + virtual void copy( TpetraVector const & x ) override; + + virtual void axpy( real64 const alpha, + TpetraVector const & x ) override; + + virtual void axpby( real64 const alpha, + TpetraVector const & x, + real64 const beta ) override; + + virtual real64 norm1() const override; + + virtual real64 norm2() const override; + + virtual real64 normInf() const override; + + virtual globalIndex globalSize() const override; + + virtual localIndex localSize() const override; + + virtual globalIndex ilower() const override; + + virtual globalIndex iupper() const override; + + virtual real64 get( globalIndex globalRow ) const override; + + virtual void get( arraySlice1d< globalIndex const > const & globalIndices, + arraySlice1d< real64 > const & values ) const override; + + virtual localIndex getLocalRowID( globalIndex const globalRow ) const override; + + virtual globalIndex getGlobalRowID( localIndex const localRow ) const override; + + virtual real64 const * extractLocalVector() const override; + + virtual real64 * extractLocalVector() override; + + virtual void extract( arrayView1d< real64 > const & localVector ) const override; + + virtual MPI_Comm getComm() const override; + + virtual void print( std::ostream & os = std::cout ) const override; + + virtual void write( string const & filename, + LAIOutputFormat const format = LAIOutputFormat::MATRIX_MARKET ) const override; + + ///@} + + /** + * @brief Alias for Tpetra map template instantiation used by this class. + */ + using Tpetra_Map = Tpetra::Map< int, globalIndex >; + + /** + * @brief Alias for specific Tpetra vector template instantiation wrapped by this class. + * + * @note This uses Tpetra's default execution/memory space. When built with CUDA support, + * this will be equal to Kokkos::Cuda, so we won't be able to create a host-only vector. + * If we want both in the same executable, we'll have to make adjustments to our LAI approach. + */ + using Tpetra_Vector = Tpetra::Vector< real64, int, globalIndex >; + + /** + * @brief Alias for specific Tpetra::MultiVector vector template instantiation wrapped by this class. + * + * This is needed for correct specification of template arguments for solver classes (e.g. Belos) + * which must be templated on MultiVector and not Vector to invoke the right explicit instantiations. + */ + using Tpetra_MultiVector = Tpetra::MultiVector< real64, int, globalIndex >; + + /** + * @brief Get the underlying Tpetra object. + * @return reference to wrapped vector + */ + Tpetra_Vector const & unwrapped() const; + + /** + * @copydoc unwrapped() const + */ + Tpetra_Vector & unwrapped(); + +private: + + /// Pointer to wrapped Tpetra object. + std::unique_ptr< Tpetra_Vector > m_vector; + +}; + +} // namespace geosx + +#endif //GEOSX_LINEARALGEBRA_INTERFACES_TPETRAVECTOR_HPP_ diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosInterface.hpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosInterface.hpp index 8d6d5b23865..bd79e1352a2 100644 --- a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosInterface.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosInterface.hpp @@ -23,6 +23,7 @@ #include "linearAlgebra/interfaces/trilinos/EpetraMatrix.hpp" #include "linearAlgebra/interfaces/trilinos/TrilinosSolver.hpp" #include "linearAlgebra/solvers/PreconditionerBase.hpp" +#include "linearAlgebra/utilities/LinearSolverParameters.hpp" #include @@ -36,7 +37,7 @@ namespace geosx struct TrilinosInterface { /** - * @brief Initializes the MPI environment for the Trilinos library + * @brief Initializes the Trilinos library * * @param[in] argc standard argc as in any C main * @param[in] argv standard argv as in any C main @@ -44,7 +45,7 @@ struct TrilinosInterface static void initialize( int & argc, char * * & argv ); /** - * @brief Finalizes the MPI environment for the Trilinos library + * @brief Finalizes the Trilinos library */ static void finalize(); diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosPreconditioner.cpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosPreconditioner.cpp index ff5b1e687b5..ef187499d7c 100644 --- a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosPreconditioner.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosPreconditioner.cpp @@ -53,33 +53,30 @@ string const & getMLCycleType( string const & value ) return optionMap.at( value ); } -string const & getMLSmootherType( string const & value ) +string const & getMLSmootherType( LinearSolverParameters::PreconditionerType const & value ) { - static std::map< string, string > const optionMap = + static std::map< LinearSolverParameters::PreconditionerType, string > const optionMap = { - { "jacobi", "Jacobi" }, - { "blockJacobi", "block Jacobi" }, - { "gaussSeidel", "Gauss-Seidel" }, - { "blockGaussSeidel", "block Gauss-Seidel" }, - { "chebyshev", "Chebyshev" }, - { "icc", "IC" }, - { "ilu", "ILU" }, - { "ilut", "ILUT" }, + { LinearSolverParameters::PreconditionerType::jacobi, "Jacobi" }, + { LinearSolverParameters::PreconditionerType::gs, "Gauss-Seidel" }, + { LinearSolverParameters::PreconditionerType::chebyshev, "Chebyshev" }, + { LinearSolverParameters::PreconditionerType::icc, "IC" }, + { LinearSolverParameters::PreconditionerType::iluk, "ILU" }, + { LinearSolverParameters::PreconditionerType::ilut, "ILUT" }, }; GEOSX_LAI_ASSERT_MSG( optionMap.count( value ) > 0, "Unsupported Trilinos/ML smoother option: " << value ); return optionMap.at( value ); } -string const & getMLCoarseType( string const & value ) +string const & getMLCoarseType( LinearSolverParameters::PreconditionerType const & value ) { - static std::map< string, string > const optionMap = + static std::map< LinearSolverParameters::PreconditionerType, string > const optionMap = { - { "jacobi", "Jacobi" }, - { "gaussSeidel", "Gauss-Seidel" }, - { "blockGaussSeidel", "block Gauss-Seidel" }, - { "chebyshev", "Chebyshev" }, - { "direct", "Amesos-KLU"}, + { LinearSolverParameters::PreconditionerType::jacobi, "Jacobi" }, + { LinearSolverParameters::PreconditionerType::gs, "Gauss-Seidel" }, + { LinearSolverParameters::PreconditionerType::chebyshev, "Chebyshev" }, + { LinearSolverParameters::PreconditionerType::direct, "Amesos-KLU"}, }; GEOSX_LAI_ASSERT_MSG( optionMap.count( value ) > 0, "Unsupported Trilinos/ML coarse solver option: " << value ); @@ -117,9 +114,11 @@ Ifpack::EPrecType getIfpackPrecondType( LinearSolverParameters::PreconditionerTy { LinearSolverParameters::PreconditionerType::ilut, Ifpack::ILUT }, { LinearSolverParameters::PreconditionerType::icc, Ifpack::IC }, { LinearSolverParameters::PreconditionerType::ict, Ifpack::ICT }, + { LinearSolverParameters::PreconditionerType::chebyshev, Ifpack::CHEBYSHEV }, { LinearSolverParameters::PreconditionerType::jacobi, Ifpack::POINT_RELAXATION }, { LinearSolverParameters::PreconditionerType::gs, Ifpack::POINT_RELAXATION }, { LinearSolverParameters::PreconditionerType::sgs, Ifpack::POINT_RELAXATION }, + { LinearSolverParameters::PreconditionerType::direct, Ifpack::AMESOS }, }; GEOSX_LAI_ASSERT_MSG( typeMap.count( type ) > 0, "Unsupported Trilinos/Ifpack preconditioner option: " << type ); @@ -146,15 +145,33 @@ CreateIfpackOperator( LinearSolverParameters const & params, Epetra_RowMatrix co const_cast< Epetra_RowMatrix * >( &matrix ), params.dd.overlap ) ); Teuchos::ParameterList list; - if( type == Ifpack::POINT_RELAXATION ) + switch( type ) { - list.set( "relaxation: type", getIfpackRelaxationType( params.preconditionerType ) ); - } - else - { - list.set( "fact: level-of-fill", params.ilu.fill ); - list.set( "fact: ilut level-of-fill", std::max( real64( params.ilu.fill ), 1.0 ) ); - list.set( "fact: ict level-of-fill", std::max( real64( params.ilu.fill ), 1.0 ) ); + case Ifpack::POINT_RELAXATION: + { + list.set( "relaxation: type", getIfpackRelaxationType( params.preconditionerType ) ); + break; + } + case Ifpack::ILU: + case Ifpack::IC: + { + list.set( "fact: level-of-fill", params.ilu.fill ); + break; + } + case Ifpack::ILUT: + { + list.set( "fact: ilut level-of-fill", std::max( real64( params.ilu.fill ), 1.0 ) ); + break; + } + case Ifpack::ICT: + { + list.set( "fact: ict level-of-fill", std::max( real64( params.ilu.fill ), 1.0 ) ); + break; + } + default: + { + break; + } } GEOSX_LAI_CHECK_ERROR( precond->SetParameters( list ) ); @@ -188,6 +205,7 @@ void TrilinosPreconditioner::compute( Matrix const & mat ) case LinearSolverParameters::PreconditionerType::ilut: case LinearSolverParameters::PreconditionerType::icc: case LinearSolverParameters::PreconditionerType::ict: + case LinearSolverParameters::PreconditionerType::direct: { m_precond = CreateIfpackOperator( m_parameters, mat.unwrapped() ); break; diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosPreconditioner.hpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosPreconditioner.hpp index ce36fee180a..d411d6e7559 100644 --- a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosPreconditioner.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosPreconditioner.hpp @@ -16,8 +16,8 @@ * @file TrilinosPreconditioner.hpp */ -#ifndef GEOSX_LINEARALGEBRA_INTERFACES_TRILINOSAMG_HPP_ -#define GEOSX_LINEARALGEBRA_INTERFACES_TRILINOSAMG_HPP_ +#ifndef GEOSX_LINEARALGEBRA_INTERFACES_TRILINOSPRECONDITIONER_HPP_ +#define GEOSX_LINEARALGEBRA_INTERFACES_TRILINOSPRECONDITIONER_HPP_ #include "linearAlgebra/solvers/PreconditionerBase.hpp" #include "linearAlgebra/interfaces/trilinos/TrilinosInterface.hpp" @@ -104,4 +104,4 @@ class TrilinosPreconditioner final : public PreconditionerBase< TrilinosInterfac } -#endif //GEOSX_LINEARALGEBRA_INTERFACES_TRILINOSAMG_HPP_ +#endif //GEOSX_LINEARALGEBRA_INTERFACES_TRILINOSPRECONDITIONER_HPP_ diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosSolver.cpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosSolver.cpp index 3113e08651e..a0c95baf059 100644 --- a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosSolver.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosSolver.cpp @@ -90,10 +90,12 @@ void TrilinosSolver::solve_direct( EpetraMatrix & mat, &rhs.unwrapped() ); // Instantiate the Amesos solver. - Amesos Factory; + Amesos factory; - // Select KLU solver (only one available as of 9/20/2018) - std::unique_ptr< Amesos_BaseSolver > solver( Factory.Create( "Klu", problem ) ); + // Select KLU solver. + // Current Trilinos release (13.0) is not compatible with our SuperLU_Dist version (6.3) + // Amesos' Umfpack interface does not work correctly with 64-bit global indices. + std::unique_ptr< Amesos_BaseSolver > solver( factory.Create( "Klu", problem ) ); // Factorize the matrix GEOSX_LAI_CHECK_ERROR( solver->SymbolicFactorization() ); diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosSolver.hpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosSolver.hpp index 8f3b89a2e54..76f441a2205 100644 --- a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosSolver.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosSolver.hpp @@ -50,13 +50,11 @@ class TrilinosSolver ~TrilinosSolver(); /** - * @brief Solve system with an iterative solver. + * @brief Solve a linear system Ax=b with an iterative solver. * @param[in,out] mat the matrix * @param[in,out] sol the solution * @param[in,out] rhs the right-hand side * @param dofManager the Degree-of-Freedom manager associated with matrix - * - * Solve Ax=b with A an EpetraMatrix, x and b EpetraVector. */ void solve( EpetraMatrix & mat, EpetraVector & sol, diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosTpetraInterface.cpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosTpetraInterface.cpp new file mode 100644 index 00000000000..7759785c4d1 --- /dev/null +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosTpetraInterface.cpp @@ -0,0 +1,45 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2019 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2019 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2019 Total, S.A + * Copyright (c) 2019- GEOSX Contributors + * All right reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file TrilinosTpetraInterface.cpp + */ + +#include "TrilinosTpetraInterface.hpp" + +#include "linearAlgebra/interfaces/trilinos/TrilinosTpetraPreconditioner.hpp" + +#include + +namespace geosx +{ + +void TrilinosTpetraInterface::initialize( int & argc, + char * * & argv ) +{ + Tpetra::initialize( &argc, &argv ); +} + +void TrilinosTpetraInterface::finalize() +{ + Tpetra::finalize(); +} + +std::unique_ptr< PreconditionerBase< TrilinosTpetraInterface > > +TrilinosTpetraInterface::createPreconditioner( LinearSolverParameters params ) +{ + return std::make_unique< TrilinosTpetraPreconditioner >( params ); +} + +} diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosTpetraInterface.hpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosTpetraInterface.hpp new file mode 100644 index 00000000000..5646e73e55a --- /dev/null +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosTpetraInterface.hpp @@ -0,0 +1,67 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2019 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2019 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2019 Total, S.A + * Copyright (c) 2019- GEOSX Contributors + * All right reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file TrilinosTpetraInterface.hpp + */ + +#ifndef GEOSX_LINEARALGEBRA_INTERFACES_TRILINOSTPETRAINTERFACE_HPP_ +#define GEOSX_LINEARALGEBRA_INTERFACES_TRILINOSTPETRAINTERFACE_HPP_ + +#include "linearAlgebra/interfaces/trilinos/TpetraVector.hpp" +#include "linearAlgebra/interfaces/trilinos/TpetraMatrix.hpp" +#include "linearAlgebra/interfaces/trilinos/TrilinosTpetraSolver.hpp" +#include "linearAlgebra/solvers/PreconditionerBase.hpp" +#include "linearAlgebra/utilities/LinearSolverParameters.hpp" + +namespace geosx +{ + +/** + * @brief This class holds aliases based on the Trilinos library. + */ +struct TrilinosTpetraInterface +{ +/** + * @brief Initializes the Trilinos library + * + * @param[in] argc standard argc as in any C main + * @param[in] argv standard argv as in any C main + */ + static void initialize( int & argc, char * * & argv ); + + /** + * @brief Finalizes the Trilinos library + */ + static void finalize(); + + /** + * @brief Create a Trilinos-based preconditioner object. + * @param params the preconditioner parameters + * @return an owning pointer to the newly created preconditioner + */ + static std::unique_ptr< PreconditionerBase< TrilinosTpetraInterface > > + createPreconditioner( LinearSolverParameters params ); + + /// Alias for TpetraMatrix + using ParallelMatrix = TpetraMatrix; + /// Alias for TpetraVector + using ParallelVector = TpetraVector; + /// Alias for TrilinosTpetraSolver + using LinearSolver = TrilinosTpetraSolver; +}; + +} + +#endif //GEOSX_LINEARALGEBRA_INTERFACES_TRILINOSTPETRAINTERFACE_HPP_ diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosTpetraPreconditioner.cpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosTpetraPreconditioner.cpp new file mode 100644 index 00000000000..2dae397f2c0 --- /dev/null +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosTpetraPreconditioner.cpp @@ -0,0 +1,260 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2019 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2019 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2019 Total, S.A + * Copyright (c) 2019- GEOSX Contributors + * All right reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file TrilinosTpetraPreconditioner.cpp + */ + +#include "TrilinosTpetraPreconditioner.hpp" + +#include +#include +#include + +namespace geosx +{ + +TrilinosTpetraPreconditioner::TrilinosTpetraPreconditioner( LinearSolverParameters params ) + : Base{}, + m_parameters( std::move( params ) ), + m_precond{} +{ } + +TrilinosTpetraPreconditioner::~TrilinosTpetraPreconditioner() = default; + +namespace +{ + +string const & getMueLuCoarseType( LinearSolverParameters::PreconditionerType const & value ) +{ + static std::map< LinearSolverParameters::PreconditionerType, string > const optionMap = + { + { LinearSolverParameters::PreconditionerType::direct, "KLU2" }, + }; + + GEOSX_LAI_ASSERT_MSG( optionMap.count( value ) > 0, "Unsupported Trilinos/MueLu coarse solver option: " << value ); + return optionMap.at( value ); +} + +string const & getMueLuVerbosity( integer const value ) +{ + static std::array< string, 5 > const optionMap = + { + "none", + "low", + "medium", + "high", + "extreme" + }; + + return optionMap[ std::max( std::min( value, integer( optionMap.size() - 1 ) ), 0 ) ]; +} + +string getIfpack2PrecondType( LinearSolverParameters::PreconditionerType const & type ) +{ + static std::map< LinearSolverParameters::PreconditionerType, string > const typeMap = + { + { LinearSolverParameters::PreconditionerType::iluk, "RILUK" }, + { LinearSolverParameters::PreconditionerType::ilut, "ILUT" }, + { LinearSolverParameters::PreconditionerType::jacobi, "RELAXATION" }, + { LinearSolverParameters::PreconditionerType::gs, "RELAXATION" }, + { LinearSolverParameters::PreconditionerType::sgs, "RELAXATION" }, + { LinearSolverParameters::PreconditionerType::direct, "AMESOS2" } + }; + + GEOSX_LAI_ASSERT_MSG( typeMap.count( type ) > 0, "Unsupported Trilinos/Ifpack2 preconditioner option: " << type ); + return typeMap.at( type ); +} + +string getIfpack2RelaxationType( LinearSolverParameters::PreconditionerType const & type ) +{ + static std::map< LinearSolverParameters::PreconditionerType, string > const typeMap = + { + { LinearSolverParameters::PreconditionerType::jacobi, "Jacobi" }, + { LinearSolverParameters::PreconditionerType::gs, "Gauss-Seidel" }, + { LinearSolverParameters::PreconditionerType::sgs, "Symmetric Gauss-Seidel" }, + }; + + GEOSX_LAI_ASSERT_MSG( typeMap.count( type ) > 0, "Unsupported Trilinos/Ifpack2 preconditioner option: " << type ); + return typeMap.at( type ); +} + +std::unique_ptr< TrilinosTpetraPreconditioner::Tpetra_Operator > +CreateMueLuOperator( LinearSolverParameters const & params, + TrilinosTpetraPreconditioner::Tpetra_RowMatrix const & matrix ) +{ + using Tpetra_Operator = TrilinosTpetraPreconditioner::Tpetra_Operator; + using Tpetra_RowMatrix = TrilinosTpetraPreconditioner::Tpetra_RowMatrix; + + Teuchos::ParameterList list; + list.set( "verbosity", getMueLuVerbosity( params.logLevel ) ); + list.set( "problem: symmetric", params.isSymmetric ); + list.set( "max levels", params.amg.maxLevels ); + list.set( "cycle type", params.amg.cycleType ); + list.set( "number of equations", params.dofsPerNode ); + list.set( "smoother: pre or post", params.amg.preOrPostSmoothing ); + list.set( "smoother: type", getIfpack2PrecondType( params.amg.smootherType ) ); + list.set( "coarse: type", getMueLuCoarseType( params.amg.coarseType ) ); + list.set( "smoother: overlap", params.dd.overlap ); + list.set( "aggregation: drop tol", params.amg.threshold ); + + Teuchos::ParameterList smootherList; + smootherList.set( "relaxation: type", getIfpack2RelaxationType( params.amg.smootherType ) ); + smootherList.set( "relaxation: sweeps", params.amg.numSweeps ); + smootherList.set( "relaxation: use l1", true ); + list.set( "smoother: params", smootherList ); + + // Const cast unfortunately required, see e.g. https://github.com/trilinos/Trilinos/issues/6390 + // TODO: upon upgrading Trilinos release, check if this has been fixed + Teuchos::RCP< Tpetra_Operator > matRCP( const_cast< Tpetra_RowMatrix * >( &matrix ), false ); + + std::unique_ptr< Tpetra_Operator > precond( MueLu::CreateTpetraPreconditioner( matRCP, list ).release().get() ); + + return precond; +} + +std::unique_ptr< TrilinosTpetraPreconditioner::Tpetra_Operator > +CreateIfpackOperator( LinearSolverParameters const & params, + TrilinosTpetraPreconditioner::Tpetra_RowMatrix const & matrix ) +{ + using Tpetra_Operator = TrilinosTpetraPreconditioner::Tpetra_Operator; + using Ifpack2_Preconditioner = Ifpack2::Preconditioner< Tpetra_Operator::scalar_type, + Tpetra_Operator::local_ordinal_type, + Tpetra_Operator::global_ordinal_type, + Tpetra_Operator::node_type >; + + Teuchos::RCP< Ifpack2_Preconditioner > precond = Ifpack2::Factory::create( "SCHWARZ", Teuchos::rcp( &matrix, false ) ); + + string const innerPrecType = getIfpack2PrecondType( params.preconditionerType ); + + Teuchos::ParameterList innerParams; + + if( innerPrecType == "RELAXATION" ) + { + innerParams.set( "relaxation: type", getIfpack2RelaxationType( params.preconditionerType ) ); + } + else if( innerPrecType == "RILUK" ) + { + innerParams.set( "fact: iluk level-of-fill", params.ilu.fill ); + innerParams.set( "fact: drop tolerance", params.ilu.threshold ); + } + else if( innerPrecType == "ILUT" ) + { + innerParams.set( "fact: ilut level-of-fill", std::fmax( params.ilu.fill, 1.0 ) ); + innerParams.set( "fact: drop tolerance", params.ilu.threshold ); + } + + Teuchos::ParameterList outerParams; + outerParams.set( "schwarz: overlap level", params.dd.overlap ); + outerParams.set( "schwarz: subdomain solver name", innerPrecType ); + outerParams.set( "schwarz: subdomain solver parameters", innerParams ); + precond->setParameters( outerParams ); + + precond->compute(); + + return std::unique_ptr< Tpetra_Operator >( precond.release().get() ); +} + +std::unique_ptr< TrilinosTpetraPreconditioner::Tpetra_Operator > +CreateIfpackDirect( LinearSolverParameters const & GEOSX_UNUSED_PARAM( params ), + TrilinosTpetraPreconditioner::Tpetra_RowMatrix const & matrix ) +{ + using Tpetra_Operator = TrilinosTpetraPreconditioner::Tpetra_Operator; + using Ifpack2_Preconditioner = Ifpack2::Preconditioner< Tpetra_Operator::scalar_type, + Tpetra_Operator::local_ordinal_type, + Tpetra_Operator::global_ordinal_type, + Tpetra_Operator::node_type >; + + Teuchos::RCP< Ifpack2_Preconditioner > precond = Ifpack2::Factory::create( "AMESOS2", Teuchos::rcp( &matrix, false ) ); + precond->compute(); + + return std::unique_ptr< Tpetra_Operator >( precond.release().get() ); +} + +} // namespace + + +void TrilinosTpetraPreconditioner::compute( Matrix const & mat ) +{ + Base::compute( mat ); + + switch( m_parameters.preconditionerType ) + { + case LinearSolverParameters::PreconditionerType::amg: + { + m_precond = CreateMueLuOperator( m_parameters, mat.unwrapped() ); + break; + } + case LinearSolverParameters::PreconditionerType::iluk: + case LinearSolverParameters::PreconditionerType::ilut: + { + m_precond = CreateIfpackOperator( m_parameters, mat.unwrapped() ); + break; + } + case LinearSolverParameters::PreconditionerType::direct: + { + m_precond = CreateIfpackDirect( m_parameters, mat.unwrapped() ); + break; + } + case LinearSolverParameters::PreconditionerType::none: + { + m_precond.reset(); + break; + } + default: + { + GEOSX_ERROR( "Unsupported preconditioner type: " << m_parameters.preconditionerType ); + } + } +} + +void TrilinosTpetraPreconditioner::apply( Vector const & src, + Vector & dst ) const +{ + GEOSX_LAI_ASSERT( ready() ); + GEOSX_LAI_ASSERT( src.ready() ); + GEOSX_LAI_ASSERT( dst.ready() ); + GEOSX_LAI_ASSERT_EQ( src.globalSize(), this->numGlobalCols() ); + GEOSX_LAI_ASSERT_EQ( dst.globalSize(), this->numGlobalRows() ); + + if( m_parameters.preconditionerType == LinearSolverParameters::PreconditionerType::none ) + { + dst.copy( src ); + } + else + { + GEOSX_LAI_ASSERT( m_precond ); + m_precond->apply( src.unwrapped(), dst.unwrapped() ); + } +} + +void TrilinosTpetraPreconditioner::clear() +{ + PreconditionerBase::clear(); + m_precond.reset(); +} + +TrilinosTpetraPreconditioner::Tpetra_Operator const & TrilinosTpetraPreconditioner::unwrapped() const +{ + GEOSX_LAI_ASSERT( ready() ); + return *m_precond; +} + +TrilinosTpetraPreconditioner::Tpetra_Operator & TrilinosTpetraPreconditioner::unwrapped() +{ + GEOSX_LAI_ASSERT( ready() ); + return *m_precond; +} + +} // namespace geosx diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosTpetraPreconditioner.hpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosTpetraPreconditioner.hpp new file mode 100644 index 00000000000..94d2aa16895 --- /dev/null +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosTpetraPreconditioner.hpp @@ -0,0 +1,109 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2019 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2019 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2019 Total, S.A + * Copyright (c) 2019- GEOSX Contributors + * All right reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file TrilinosTpetraPreconditioner.hpp + */ + +#ifndef GEOSX_LINEARALGEBRA_INTERFACES_TRILINOSTPETRAPRECONDITIONER_HPP_ +#define GEOSX_LINEARALGEBRA_INTERFACES_TRILINOSTPETRAPRECONDITIONER_HPP_ + +#include "linearAlgebra/solvers/PreconditionerBase.hpp" +#include "linearAlgebra/interfaces/trilinos/TrilinosTpetraInterface.hpp" +#include "linearAlgebra/utilities/LinearSolverParameters.hpp" + +#include +#include + +#include + +namespace geosx +{ + +/** + * @brief Wrapper around Trilinos Tpetra-based preconditioners. + */ +class TrilinosTpetraPreconditioner final : public PreconditionerBase< TrilinosTpetraInterface > +{ +public: + + /// Alias for base type + using Base = PreconditionerBase< TrilinosTpetraInterface >; + + /// Alias for vector type + using Vector = typename Base::Vector; + + /// Alias for matrix type + using Matrix = typename Base::Matrix; + + /// Allow for partial overload of Base::compute() + using Base::compute; + + /** + * @brief Constructor. + * @param params preconditioner parameters + */ + explicit TrilinosTpetraPreconditioner( LinearSolverParameters params ); + + /** + * @brief Destructor. + */ + virtual ~TrilinosTpetraPreconditioner() override; + + /** + * @brief Compute the preconditioner from a matrix. + * @param mat the matrix to precondition. + */ + virtual void compute( Matrix const & mat ) override; + + /** + * @brief Apply operator to a vector + * @param src Input vector (x). + * @param dst Output vector (b). + * + * @warning @p src and @p dst cannot alias the same vector. + */ + virtual void apply( Vector const & src, Vector & dst ) const override; + + virtual void clear() override; + + /// Instantiation of Tpetra::Operator template held by this class + using Tpetra_Operator = Tpetra::Operator< real64, int, globalIndex >; + + /// Instantiation of Tpetra::RowMatrix template used by this class + using Tpetra_RowMatrix = Tpetra::RowMatrix< real64, int, globalIndex >; + + /** + * @brief Access the underlying Epetra preconditioning operator. + * @return reference to the Epetra operator + */ + Tpetra_Operator const & unwrapped() const; + + /** + * @copydoc unwrapped() const + */ + Tpetra_Operator & unwrapped(); + +private: + + /// Parameters for all preconditioners + LinearSolverParameters m_parameters; + + /// Pointer to the Trilinos implementation + std::unique_ptr< Tpetra_Operator > m_precond; +}; + +} // namespace geosx + +#endif //GEOSX_LINEARALGEBRA_INTERFACES_TRILINOSTPETRAPRECONDITIONER_HPP_ diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosTpetraSolver.cpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosTpetraSolver.cpp new file mode 100644 index 00000000000..f99ff61f36f --- /dev/null +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosTpetraSolver.cpp @@ -0,0 +1,291 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2019 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2019 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2019 Total, S.A + * Copyright (c) 2019- GEOSX Contributors + * All right reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file TrilinosTpetraSolver.cpp + */ + +#include "TrilinosTpetraSolver.hpp" + +#include "common/Stopwatch.hpp" +#include "linearAlgebra/DofManager.hpp" +#include "linearAlgebra/interfaces/trilinos/TpetraMatrix.hpp" +#include "linearAlgebra/interfaces/trilinos/TpetraVector.hpp" +#include "linearAlgebra/interfaces/trilinos/TrilinosTpetraPreconditioner.hpp" +#include "linearAlgebra/utilities/LAIHelperFunctions.hpp" + +#include +#include +#include + +namespace geosx +{ + +TrilinosTpetraSolver::TrilinosTpetraSolver( LinearSolverParameters parameters ): + m_parameters( std::move( parameters ) ) +{} + +TrilinosTpetraSolver::~TrilinosTpetraSolver() = default; + +void TrilinosTpetraSolver::solve( TpetraMatrix & mat, + TpetraVector & sol, + TpetraVector & rhs, + DofManager const * const GEOSX_UNUSED_PARAM( dofManager ) ) +{ + GEOSX_LAI_ASSERT( mat.ready() ); + GEOSX_LAI_ASSERT( sol.ready() ); + GEOSX_LAI_ASSERT( rhs.ready() ); + + if( m_parameters.solverType == LinearSolverParameters::SolverType::direct ) + { + solve_direct( mat, sol, rhs ); + } + else + { + solve_krylov( mat, sol, rhs ); + } +} + +namespace +{ + +string getAmesos2SolverName( bool const parallel ) +{ + // Define a list of direct solvers to try ordered by priority + static std::array< string, 1 > parallelSolvers = + { + "superlu_dist", + }; + + static std::array< string, 2 > serialSolvers = + { + "umfpack", + "Klu2" + }; + + // First try parallel solvers, if requested (which is normally the default) + if( parallel ) + { + for( string const & solver : parallelSolvers ) + { + if( Amesos2::query( solver ) ) + { + return solver; + } + } + } + + // Fall back to serial if parallel not available or user doesn't want it + for( string const & solver : serialSolvers ) + { + if( Amesos2::query( solver ) ) + { + return solver; + } + } + GEOSX_ERROR( "No direct solver available in Trilinos/Amesos2" ); + return ""; +} + +string getAmesos2SuperLUDistColPerm( LinearSolverParameters::Direct::ColPerm const value ) +{ + static std::map< LinearSolverParameters::Direct::ColPerm, string > const optionMap = + { + { LinearSolverParameters::Direct::ColPerm::none, "NATURAL" }, + { LinearSolverParameters::Direct::ColPerm::parmetis, "PARMETIS" }, + }; + + // Since not all valid options can be accepted by Amesos2, return default value instead of error. + string result; + if( optionMap.count( value ) > 0 ) + { + result = optionMap.at( value ); + } + else + { + LinearSolverParameters::Direct::ColPerm const defaultValue = LinearSolverParameters::Direct::ColPerm::none; + GEOSX_WARNING( "Unsupported Trilinos/Amesos2/SuperLU_Dist column permutation option: " << value << ". " + "Default value '" << defaultValue << "' will be used." ); + result = optionMap.at( defaultValue ); + } + return result; +} + +} // namespace + +void TrilinosTpetraSolver::solve_direct( TpetraMatrix & mat, + TpetraVector & sol, + TpetraVector & rhs ) +{ + // Time setup and solve + Stopwatch watch; + + // Create Epetra linear problem and instantiate solver. + using Tpetra_Matrix = TpetraMatrix::Tpetra_CrsMatrix; + using Tpetra_MultiVector = TpetraVector::Tpetra_MultiVector; + using Amesos2_Solver = Amesos2::Solver< Tpetra_Matrix, Tpetra_MultiVector >; + + // Instantiate the solver + // Note: we explicitly provide template arguments to create() because + // Tpetra::Vector is not actually recognized by some Amesos2 internals + string const solverName = getAmesos2SolverName( m_parameters.direct.parallel ); + Teuchos::RCP< Amesos2_Solver > solver = + Amesos2::create< Tpetra_Matrix, Tpetra_MultiVector >( solverName.c_str(), + Teuchos::rcp( &mat.unwrapped(), false ), + Teuchos::rcp( &sol.unwrapped(), false ), + Teuchos::rcp( &rhs.unwrapped(), false ) ); + + // Set additional solver options + Teuchos::ParameterList params( "Amesos2" ); + if( solverName == "superlu_dist" ) + { + // Only these two options are supported by Amesos2_Superludist + params.set( "ColPerm", getAmesos2SuperLUDistColPerm( m_parameters.direct.colPerm ) ); + params.set( "ReplaceTinyPivot", bool( m_parameters.direct.replaceTinyPivot ) ); + } + solver->setParameters( Teuchos::rcp( ¶ms, false ) ); + + // Factorize the matrix + solver->symbolicFactorization(); + solver->numericFactorization(); + m_result.setupTime = watch.elapsedTime(); + + // Solve the system + watch.zero(); + solver->solve(); + m_result.solveTime = watch.elapsedTime(); + + // Basic output + if( m_parameters.logLevel > 0 ) + { + solver->printTiming( *Teuchos::getFancyOStream( Teuchos::rcp( &std::cout, false ) ) ); + } + + m_result.status = LinearSolverResult::Status::Success; + m_result.numIterations = 1; + m_result.residualReduction = NumericTraits< real64 >::eps; +} + +namespace +{ + +string getBelosSolverName( LinearSolverParameters::SolverType const & type ) +{ + static std::map< LinearSolverParameters::SolverType, string > typeMap = + { + { LinearSolverParameters::SolverType::cg, "CG" }, + { LinearSolverParameters::SolverType::bicgstab, "BICGSTAB" }, + { LinearSolverParameters::SolverType::gmres, "GMRES" }, + { LinearSolverParameters::SolverType::fgmres, "FLEXIBLE GMRES" } + }; + + GEOSX_LAI_ASSERT_MSG( typeMap.count( type ) > 0, "Unsupported Trilinos/Belos solver option: " << type ); + return typeMap.at( type ); +} + +Belos::MsgType getBelosVerbosity( integer const value ) +{ + static std::array< Belos::MsgType, 5 > const optionMap = + { + Belos::MsgType( Belos::Errors ), + Belos::MsgType( Belos::Errors | Belos::Warnings | Belos::FinalSummary ), + Belos::MsgType( Belos::Errors | Belos::Warnings | Belos::FinalSummary | Belos::IterationDetails ), + Belos::MsgType( Belos::Errors | Belos::Warnings | Belos::FinalSummary | Belos::IterationDetails | Belos::TimingDetails ), + Belos::MsgType( Belos::Errors | Belos::Warnings | Belos::FinalSummary | Belos::IterationDetails | Belos::TimingDetails | Belos::Debug ) + }; + + return optionMap[ std::max( std::min( value, integer( optionMap.size() - 1 ) ), 0 ) ]; +} + +integer getBelosOutputFrequency( integer const value ) +{ + static std::array< integer, 5 > const optionMap = + { + -1, + 1, + 1, + 1, + 1 + }; + + return optionMap[ std::max( std::min( value, integer( optionMap.size() - 1 ) ), 0 ) ]; +} + +} // namespace + +void TrilinosTpetraSolver::solve_krylov( TpetraMatrix & mat, + TpetraVector & sol, + TpetraVector & rhs ) +{ + // Time setup and solve + Stopwatch watch; + + using Tpetra_MultiVector = TpetraVector::Tpetra_MultiVector; + using Tpetra_Operator = TrilinosTpetraPreconditioner::Tpetra_Operator; + using Belos_LinearProblem = Belos::LinearProblem< real64, Tpetra_MultiVector, Tpetra_Operator >; + using Belos_Factory = Belos::SolverFactory< real64, Tpetra_MultiVector, Tpetra_Operator >; + using Belos_Solver = Belos::SolverManager< real64, Tpetra_MultiVector, Tpetra_Operator >; + + // Deal with separate component approximation + TpetraMatrix separateComponentMatrix; + if( m_parameters.amg.separateComponents ) + { + LAIHelperFunctions::SeparateComponentFilter( mat, separateComponentMatrix, m_parameters.dofsPerNode ); + } + TpetraMatrix & precondMat = m_parameters.amg.separateComponents ? separateComponentMatrix : mat; + + // Create and compute preconditioner + TrilinosTpetraPreconditioner precond( m_parameters ); + precond.compute( precondMat ); + + // Create a linear problem. + Teuchos::RCP< Belos_LinearProblem > problem = + Teuchos::rcp( new Belos_LinearProblem( Teuchos::rcp( &mat.unwrapped(), false ), + Teuchos::rcp( &sol.unwrapped(), false ), + Teuchos::rcp( &rhs.unwrapped(), false ) ) ); + + // Set as right preconditioner. + // TODO: do we need an option for left/right? + problem->setRightPrec( Teuchos::rcp( &precond.unwrapped(), false ) ); + problem->setProblem(); + + Teuchos::ParameterList list; + list.set( "Implicit Residual Scaling", "Norm of RHS" ); + list.set( "Explicit Residual Scaling", "Norm of RHS" ); + list.set( "Maximum Iterations", m_parameters.krylov.maxIterations ); + list.set( "Num Blocks", m_parameters.krylov.maxRestart ); + list.set( "Maximum Restarts", m_parameters.krylov.maxIterations / m_parameters.krylov.maxRestart - 1 ); + list.set( "Orthogonalization", "DGKS" ); // this one seems more robust + list.set( "Convergence Tolerance", m_parameters.krylov.relTolerance ); + list.set( "Verbosity", getBelosVerbosity( m_parameters.logLevel ) ); + list.set( "Output Frequency", getBelosOutputFrequency( m_parameters.logLevel ) ); + + Belos_Factory factory; + Teuchos::RCP< Belos_Solver > solver = factory.create( getBelosSolverName( m_parameters.solverType ), + Teuchos::rcp( &list, false ) ); + solver->setProblem( problem ); + + m_result.setupTime = watch.elapsedTime(); + + // Do solve + watch.zero(); + Belos::ReturnType const result = solver->solve(); + m_result.solveTime = watch.elapsedTime(); + + m_result.status = ( result == Belos::Converged ) ? LinearSolverResult::Status::Success : LinearSolverResult::Status::NotConverged; + m_result.numIterations = solver->getNumIters(); + m_result.residualReduction = solver->achievedTol(); +} + +} // end geosx namespace diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosTpetraSolver.hpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosTpetraSolver.hpp new file mode 100644 index 00000000000..00da7adf919 --- /dev/null +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosTpetraSolver.hpp @@ -0,0 +1,88 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2019 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2019 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2019 Total, S.A + * Copyright (c) 2019- GEOSX Contributors + * All right reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file TrilinosTpetraSolver.hpp + */ + +#ifndef GEOSX_LINEARALGEBRA_INTERFACES_TRILINOSTPETRASOLVER_HPP_ +#define GEOSX_LINEARALGEBRA_INTERFACES_TRILINOSTPETRASOLVER_HPP_ + +#include "linearAlgebra/utilities/LinearSolverParameters.hpp" +#include "linearAlgebra/utilities/LinearSolverResult.hpp" + +namespace geosx +{ + +class DofManager; +class TpetraVector; +class TpetraMatrix; + +/** + * @brief Wrapper for Trilinos/Tpetra-based direct and interative linear solvers. + */ +class TrilinosTpetraSolver +{ +public: + + /** + * @brief Solver constructor, with parameter list reference. + * @param[in] parameters structure containing linear solver parameters + */ + TrilinosTpetraSolver( LinearSolverParameters parameters ); + + /** + * @brief Virtual destructor. + * + */ + ~TrilinosTpetraSolver(); + + /** + * @brief Solve a linear system Ax=b with an iterative solver. + * @param[in,out] mat the matrix + * @param[in,out] sol the solution + * @param[in,out] rhs the right-hand side + * @param dofManager the Degree-of-Freedom manager associated with matrix + */ + void solve( TpetraMatrix & mat, + TpetraVector & sol, + TpetraVector & rhs, + DofManager const * const dofManager = nullptr ); + + /** + * @brief Get the result of previous solve. + * @return struct with last solve stats + */ + LinearSolverResult const & result() + { + return m_result; + } + +private: + + LinearSolverParameters m_parameters; + LinearSolverResult m_result; + + void solve_direct( TpetraMatrix & mat, + TpetraVector & sol, + TpetraVector & rhs ); + + void solve_krylov( TpetraMatrix & mat, + TpetraVector & sol, + TpetraVector & rhs ); +}; + +} // namespace geosx + +#endif //GEOSX_LINEARALGEBRA_INTERFACES_TRILINOSTPETRASOLVER_HPP_ diff --git a/src/coreComponents/linearAlgebra/solvers/BiCGSTABsolver.cpp b/src/coreComponents/linearAlgebra/solvers/BiCGSTABsolver.cpp index 247e72fc07e..4cf2712ce98 100644 --- a/src/coreComponents/linearAlgebra/solvers/BiCGSTABsolver.cpp +++ b/src/coreComponents/linearAlgebra/solvers/BiCGSTABsolver.cpp @@ -160,6 +160,8 @@ void BiCGSTABsolver< VECTOR >::solve( Vector const & b, #ifdef GEOSX_USE_TRILINOS template class BiCGSTABsolver< TrilinosInterface::ParallelVector >; template class BiCGSTABsolver< BlockVectorView< TrilinosInterface::ParallelVector > >; +template class BiCGSTABsolver< TrilinosTpetraInterface::ParallelVector >; +template class BiCGSTABsolver< BlockVectorView< TrilinosTpetraInterface::ParallelVector > >; #endif #ifdef GEOSX_USE_HYPRE diff --git a/src/coreComponents/linearAlgebra/solvers/BlockPreconditioner.cpp b/src/coreComponents/linearAlgebra/solvers/BlockPreconditioner.cpp index 7d32c11170f..7f17980afc7 100644 --- a/src/coreComponents/linearAlgebra/solvers/BlockPreconditioner.cpp +++ b/src/coreComponents/linearAlgebra/solvers/BlockPreconditioner.cpp @@ -86,7 +86,7 @@ void BlockPreconditioner< LAI >::applyBlockScaling() { if( m_scalingOption == BlockScalingOption::FrobeniusNorm ) { - real64 norms[2] = { m_matBlocks( 0, 0 ).normFrobenius(), m_matBlocks( 1, 1 ).normFrobenius() }; + real64 const norms[2] = { m_matBlocks( 0, 0 ).normFrobenius(), m_matBlocks( 1, 1 ).normFrobenius() }; m_scaling[0] = std::min( norms[1] / norms[0], 1.0 ); m_scaling[1] = std::min( norms[0] / norms[1], 1.0 ); } @@ -241,6 +241,7 @@ void BlockPreconditioner< LAI >::clear() // ----------------------- #ifdef GEOSX_USE_TRILINOS template class BlockPreconditioner< TrilinosInterface >; +template class BlockPreconditioner< TrilinosTpetraInterface >; #endif #ifdef GEOSX_USE_HYPRE diff --git a/src/coreComponents/linearAlgebra/solvers/CGsolver.cpp b/src/coreComponents/linearAlgebra/solvers/CGsolver.cpp index 4f8e5092a37..a82c0e342a6 100644 --- a/src/coreComponents/linearAlgebra/solvers/CGsolver.cpp +++ b/src/coreComponents/linearAlgebra/solvers/CGsolver.cpp @@ -150,6 +150,8 @@ void CGsolver< VECTOR >::solve( Vector const & b, Vector & x ) const #ifdef GEOSX_USE_TRILINOS template class CGsolver< TrilinosInterface::ParallelVector >; template class CGsolver< BlockVectorView< TrilinosInterface::ParallelVector > >; +template class CGsolver< TrilinosTpetraInterface::ParallelVector >; +template class CGsolver< BlockVectorView< TrilinosTpetraInterface::ParallelVector > >; #endif #ifdef GEOSX_USE_HYPRE diff --git a/src/coreComponents/linearAlgebra/solvers/GMRESsolver.cpp b/src/coreComponents/linearAlgebra/solvers/GMRESsolver.cpp index 1863e5c91fa..40f75014ed0 100644 --- a/src/coreComponents/linearAlgebra/solvers/GMRESsolver.cpp +++ b/src/coreComponents/linearAlgebra/solvers/GMRESsolver.cpp @@ -209,6 +209,8 @@ void GMRESsolver< VECTOR >::solve( Vector const & b, #ifdef GEOSX_USE_TRILINOS template class GMRESsolver< TrilinosInterface::ParallelVector >; template class GMRESsolver< BlockVectorView< TrilinosInterface::ParallelVector > >; +template class GMRESsolver< TrilinosTpetraInterface::ParallelVector >; +template class GMRESsolver< BlockVectorView< TrilinosTpetraInterface::ParallelVector > >; #endif #ifdef GEOSX_USE_HYPRE diff --git a/src/coreComponents/linearAlgebra/solvers/KrylovSolver.cpp b/src/coreComponents/linearAlgebra/solvers/KrylovSolver.cpp index dbd60d4ddde..002c92f7cda 100644 --- a/src/coreComponents/linearAlgebra/solvers/KrylovSolver.cpp +++ b/src/coreComponents/linearAlgebra/solvers/KrylovSolver.cpp @@ -91,6 +91,8 @@ KrylovSolver< VECTOR >::Create( LinearSolverParameters const & parameters, #ifdef GEOSX_USE_TRILINOS template class KrylovSolver< TrilinosInterface::ParallelVector >; template class KrylovSolver< BlockVectorView< TrilinosInterface::ParallelVector > >; +template class KrylovSolver< TrilinosTpetraInterface::ParallelVector >; +template class KrylovSolver< BlockVectorView< TrilinosTpetraInterface::ParallelVector > >; #endif #ifdef GEOSX_USE_HYPRE diff --git a/src/coreComponents/linearAlgebra/solvers/SeparateComponentPreconditioner.cpp b/src/coreComponents/linearAlgebra/solvers/SeparateComponentPreconditioner.cpp index 734b741b2c2..ca89946cab8 100644 --- a/src/coreComponents/linearAlgebra/solvers/SeparateComponentPreconditioner.cpp +++ b/src/coreComponents/linearAlgebra/solvers/SeparateComponentPreconditioner.cpp @@ -70,6 +70,7 @@ void SeparateComponentPreconditioner< LAI >::clear() // ----------------------- #ifdef GEOSX_USE_TRILINOS template class SeparateComponentPreconditioner< TrilinosInterface >; +template class SeparateComponentPreconditioner< TrilinosTpetraInterface >; #endif #ifdef GEOSX_USE_HYPRE diff --git a/src/coreComponents/linearAlgebra/unitTests/CMakeLists.txt b/src/coreComponents/linearAlgebra/unitTests/CMakeLists.txt index 4ce8da4720d..92a168a5a7c 100644 --- a/src/coreComponents/linearAlgebra/unitTests/CMakeLists.txt +++ b/src/coreComponents/linearAlgebra/unitTests/CMakeLists.txt @@ -2,11 +2,13 @@ # Specify list of tests # set( LAI_tests - testLAOperations.cpp + testMatrices.cpp testArrayLAOperations.cpp + testExternalSolvers.cpp testKrylovSolvers.cpp testDofManager.cpp - testLAIHelperFunctions.cpp) + testLAIHelperFunctions.cpp + testVectors.cpp ) set( nranks 2 ) diff --git a/src/coreComponents/linearAlgebra/unitTests/testDofManager.cpp b/src/coreComponents/linearAlgebra/unitTests/testDofManager.cpp index 6dc1a52bfd0..4e32e7688d1 100644 --- a/src/coreComponents/linearAlgebra/unitTests/testDofManager.cpp +++ b/src/coreComponents/linearAlgebra/unitTests/testDofManager.cpp @@ -732,6 +732,7 @@ REGISTER_TYPED_TEST_SUITE_P( DofManagerSparsityTest, #ifdef GEOSX_USE_TRILINOS INSTANTIATE_TYPED_TEST_SUITE_P( Trilinos, DofManagerSparsityTest, TrilinosInterface, ); +INSTANTIATE_TYPED_TEST_SUITE_P( TrilinosTpetra, DofManagerSparsityTest, TrilinosTpetraInterface, ); #endif #ifdef GEOSX_USE_HYPRE diff --git a/src/coreComponents/linearAlgebra/unitTests/testExternalSolvers.cpp b/src/coreComponents/linearAlgebra/unitTests/testExternalSolvers.cpp new file mode 100644 index 00000000000..a5ac9a86c45 --- /dev/null +++ b/src/coreComponents/linearAlgebra/unitTests/testExternalSolvers.cpp @@ -0,0 +1,250 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2019 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2019 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2019 Total, S.A + * Copyright (c) 2019- GEOSX Contributors + * All right reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file testExternalSolvers.cpp + */ + +#include + +#include "common/DataTypes.hpp" +#include "linearAlgebra/interfaces/InterfaceTypes.hpp" +#include "linearAlgebra/utilities/LinearSolverParameters.hpp" +#include "managers/initialization.hpp" + +#include "testLinearAlgebraUtils.hpp" + +using namespace geosx; + +static real64 constexpr machinePrecision = 20.0 * std::numeric_limits< real64 >::epsilon(); + +/////////////////////////////////////////////////////////////////////////////////////// + +LinearSolverParameters params_Direct() +{ + LinearSolverParameters parameters; + parameters.solverType = LinearSolverParameters::SolverType::direct; + parameters.krylov.relTolerance = machinePrecision; + return parameters; +} + +LinearSolverParameters params_GMRES_ILU() +{ + LinearSolverParameters parameters; + parameters.krylov.relTolerance = 1e-8; + parameters.krylov.maxIterations = 300; + parameters.solverType = LinearSolverParameters::SolverType::gmres; + parameters.preconditionerType = LinearSolverParameters::PreconditionerType::iluk; + parameters.ilu.fill = 1; + return parameters; +} + +LinearSolverParameters params_GMRES_AMG() +{ + LinearSolverParameters parameters; + parameters.krylov.relTolerance = 1e-8; + parameters.krylov.maxIterations = 300; + parameters.solverType = LinearSolverParameters::SolverType::gmres; + parameters.preconditionerType = LinearSolverParameters::PreconditionerType::amg; + parameters.amg.smootherType = LinearSolverParameters::PreconditionerType::gs; + parameters.amg.coarseType = LinearSolverParameters::PreconditionerType::direct; + return parameters; +} + +LinearSolverParameters params_CG_AMG() +{ + LinearSolverParameters parameters; + parameters.krylov.relTolerance = 1e-8; + parameters.krylov.maxIterations = 300; + parameters.solverType = LinearSolverParameters::SolverType::cg; + parameters.isSymmetric = true; + parameters.preconditionerType = LinearSolverParameters::PreconditionerType::amg; + parameters.amg.smootherType = LinearSolverParameters::PreconditionerType::gs; + parameters.amg.coarseType = LinearSolverParameters::PreconditionerType::direct; + return parameters; +} + +template< typename LAI > +class SolverTestBase : public ::testing::Test +{ +public: + + using Matrix = typename LAI::ParallelMatrix; + using Vector = typename LAI::ParallelVector; + using Solver = typename LAI::LinearSolver; + +protected: + + Matrix matrix; + real64 cond_est = 1.0; + + void test( LinearSolverParameters const & params ) + { + // Create a random "true" solution vector + Vector sol_true; + sol_true.createWithLocalSize( matrix.numLocalCols(), matrix.getComm() ); + sol_true.rand(); + + // Create and compute the right-hand side vector + Vector rhs; + rhs.createWithLocalSize( matrix.numLocalRows(), matrix.getComm() ); + matrix.apply( sol_true, rhs ); + + // Create and zero out the computed solution vector + Vector sol_comp; + sol_comp.createWithLocalSize( sol_true.localSize(), sol_true.getComm() ); + sol_comp.zero(); + + // Create the solver and solve the system + Solver solver( params ); + solver.solve( matrix, sol_comp, rhs ); + EXPECT_TRUE( solver.result().success() ); + + // Check that solution is within epsilon of true + sol_comp.axpy( -1.0, sol_true ); + real64 const relTol = cond_est * params.krylov.relTolerance; + EXPECT_LT( sol_comp.norm2() / sol_true.norm2(), relTol ); + } +}; + +/////////////////////////////////////////////////////////////////////////////////////// + +template< typename LAI > +class SolverTestLaplace2D : public SolverTestBase< LAI > +{ +public: + + using Base = SolverTestBase< LAI >; + using Matrix = typename Base::Matrix; + using Vector = typename Base::Vector; + +protected: + + void SetUp() override + { + globalIndex constexpr n = 100; + compute2DLaplaceOperator( MPI_COMM_GEOSX, n, this->matrix ); + + // Condition number for the Laplacian matrix estimate: 4 * n^2 / pi^2 + this->cond_est = 4.0 * n * n / std::pow( M_PI, 2 ); + } +}; + +TYPED_TEST_SUITE_P( SolverTestLaplace2D ); + +TYPED_TEST_P( SolverTestLaplace2D, Direct ) +{ + this->test( params_Direct() ); +} + +TYPED_TEST_P( SolverTestLaplace2D, GMRES_ILU ) +{ + this->test( params_GMRES_ILU() ); +} + +TYPED_TEST_P( SolverTestLaplace2D, CG_AMG ) +{ + this->test( params_CG_AMG() ); +} + +REGISTER_TYPED_TEST_SUITE_P( SolverTestLaplace2D, + Direct, + GMRES_ILU, + CG_AMG ); + +#ifdef GEOSX_USE_TRILINOS +INSTANTIATE_TYPED_TEST_SUITE_P( Trilinos, SolverTestLaplace2D, TrilinosInterface, ); +INSTANTIATE_TYPED_TEST_SUITE_P( TrilinosTpetra, SolverTestLaplace2D, TrilinosTpetraInterface, ); +#endif + +#ifdef GEOSX_USE_HYPRE +INSTANTIATE_TYPED_TEST_SUITE_P( Hypre, SolverTestLaplace2D, HypreInterface, ); +#endif + +#ifdef GEOSX_USE_PETSC +INSTANTIATE_TYPED_TEST_SUITE_P( Petsc, SolverTestLaplace2D, PetscInterface, ); +#endif + +/////////////////////////////////////////////////////////////////////////////////////// + +template< typename LAI > +class SolverTestElasticity2D : public SolverTestBase< LAI > +{ +public: + + using Base = SolverTestBase< LAI >; + using Matrix = typename Base::Matrix; + using Vector = typename Base::Vector; + +protected: + + void SetUp() override + { + globalIndex constexpr n = 100; + compute2DElasticityOperator( MPI_COMM_GEOSX, 1.0, 1.0, n, n, 10000., 0.2, this->matrix ); + + // Impose Dirichlet boundary conditions: fix domain bottom (first 2*(nCellsX + 1) rows of matrix) + this->matrix.open(); + for( globalIndex iRow = 0; iRow < 2 * (n + 1); ++iRow ) + { + if( this->matrix.getLocalRowID( iRow ) >= 0 ) + { + this->matrix.clearRow( iRow, true ); + } + } + this->matrix.close(); + this->cond_est = 1e4; // not a true condition number estimate, but enough to pass tests + } +}; + +TYPED_TEST_SUITE_P( SolverTestElasticity2D ); + +TYPED_TEST_P( SolverTestElasticity2D, Direct ) +{ + this->test( params_Direct() ); +} + +TYPED_TEST_P( SolverTestElasticity2D, GMRES_AMG ) +{ + LinearSolverParameters params = params_GMRES_AMG(); + params.amg.separateComponents = true; + params.dofsPerNode = 2; + this->test( params ); +} + +REGISTER_TYPED_TEST_SUITE_P( SolverTestElasticity2D, + Direct, + GMRES_AMG ); + +#ifdef GEOSX_USE_TRILINOS +INSTANTIATE_TYPED_TEST_SUITE_P( Trilinos, SolverTestElasticity2D, TrilinosInterface, ); +INSTANTIATE_TYPED_TEST_SUITE_P( TrilinosTpetra, SolverTestElasticity2D, TrilinosTpetraInterface, ); +#endif + +#ifdef GEOSX_USE_HYPRE +INSTANTIATE_TYPED_TEST_SUITE_P( Hypre, SolverTestElasticity2D, HypreInterface, ); +#endif + +#ifdef GEOSX_USE_PETSC +INSTANTIATE_TYPED_TEST_SUITE_P( Petsc, SolverTestElasticity2D, PetscInterface, ); +#endif + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + geosx::basicSetup( argc, argv ); + int const result = RUN_ALL_TESTS(); + geosx::basicCleanup(); + return result; +} diff --git a/src/coreComponents/linearAlgebra/unitTests/testKrylovSolvers.cpp b/src/coreComponents/linearAlgebra/unitTests/testKrylovSolvers.cpp index ce7a7d5311a..8165b3fa20c 100644 --- a/src/coreComponents/linearAlgebra/unitTests/testKrylovSolvers.cpp +++ b/src/coreComponents/linearAlgebra/unitTests/testKrylovSolvers.cpp @@ -252,6 +252,7 @@ REGISTER_TYPED_TEST_SUITE_P( KrylovSolverBlockTest, #ifdef GEOSX_USE_TRILINOS INSTANTIATE_TYPED_TEST_SUITE_P( Trilinos, KrylovSolverBlockTest, TrilinosInterface, ); +INSTANTIATE_TYPED_TEST_SUITE_P( TrilinosTpetra, KrylovSolverBlockTest, TrilinosTpetraInterface, ); #endif #ifdef GEOSX_USE_HYPRE diff --git a/src/coreComponents/linearAlgebra/unitTests/testLAOperations.cpp b/src/coreComponents/linearAlgebra/unitTests/testLAOperations.cpp deleted file mode 100644 index fccadaf2642..00000000000 --- a/src/coreComponents/linearAlgebra/unitTests/testLAOperations.cpp +++ /dev/null @@ -1,799 +0,0 @@ -/* - * ------------------------------------------------------------------------------------------------------------ - * 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. - * ------------------------------------------------------------------------------------------------------------ - */ - -/** - * @file testLAOperations.cpp - */ - -#include - -#include "common/DataTypes.hpp" -#include "linearAlgebra/interfaces/InterfaceTypes.hpp" -#include "linearAlgebra/utilities/LinearSolverParameters.hpp" -#include "managers/initialization.hpp" -#include "testLinearAlgebraUtils.hpp" - -using namespace geosx; - -static real64 constexpr machinePrecision = 20.0 * std::numeric_limits< real64 >::epsilon(); - -template< typename LAI > -class LAOperationsTest : public ::testing::Test -{}; - -TYPED_TEST_SUITE_P( LAOperationsTest ); - -TYPED_TEST_P( LAOperationsTest, VectorFunctions ) -{ - // Define alias - using Vector = typename TypeParam::ParallelVector; - - // Get the MPI rank - int const rank = MpiWrapper::Comm_rank( MPI_COMM_GEOSX ); - int const numRanks = MpiWrapper::Comm_size( MPI_COMM_GEOSX ); - - Vector x; - localIndex const localSize = 3; - globalIndex const globalSize = localSize * numRanks; - globalIndex const offset = rank * localSize; - - // Testing createWithLocalSize - x.createWithLocalSize( localSize, MPI_COMM_GEOSX ); - EXPECT_EQ( x.localSize(), localSize ); - EXPECT_EQ( x.globalSize(), globalSize ); - - // Testing iupper/ilower - EXPECT_EQ( x.ilower(), offset ); - EXPECT_EQ( x.iupper(), offset + localSize ); - - // Testing setting/getting values locally - x.open(); - for( globalIndex i = x.ilower(); i < x.iupper(); ++i ) - { - x.set( i, 2 * i ); - } - x.close(); - for( globalIndex i = x.ilower(); i < x.iupper(); ++i ) - { - EXPECT_DOUBLE_EQ( x.get( i ), 2 * i ); - } - - // Testing createWithGlobalSize - x.createWithGlobalSize( globalSize, MPI_COMM_GEOSX ); - EXPECT_EQ( x.localSize(), localSize ); - EXPECT_EQ( x.globalSize(), globalSize ); - - // Testing adding global values on rank 0 and getting locally - x.open(); - if( rank == 0 ) - { - for( globalIndex i = 0; i < x.globalSize(); ++i ) - { - x.add( i, 2 * i ); - } - } - x.close(); - for( globalIndex i = x.ilower(); i < x.iupper(); ++i ) - { - EXPECT_EQ( x.get( i ), 2 * i ); - } - - // Testing getLocalRowID - for( globalIndex i = x.ilower(); i < x.iupper(); ++i ) - { - EXPECT_EQ( x.getLocalRowID( i ), i % localSize ); - } - - // Testing create with array1d - array1d< real64 > localVals( localSize ); - for( localIndex i = 0; i < localSize; ++i ) - { - localVals[i] = real64( i + rank * localSize ); - } - - Vector v; - v.create( localVals, MPI_COMM_GEOSX ); - for( globalIndex i = v.ilower(); i < v.iupper(); ++i ) - { - EXPECT_EQ( v.get( i ), localVals[v.getLocalRowID( i )] ); - } - - // Testing copy constructor, assignment, get element - Vector y( x ); - Vector z; - z = x; - for( globalIndex i = x.ilower(); i < x.iupper(); ++i ) - { - EXPECT_EQ( x.get( i ), y.get( i ) ); - EXPECT_EQ( x.get( i ), z.get( i ) ); - } - - // Testing zero - z.zero(); - for( globalIndex i = y.ilower(); i < y.iupper(); ++i ) - { - EXPECT_EQ( z.get( i ), 0 ); - } - - // Testing copy - z.copy( x ); - for( globalIndex i = y.ilower(); i < y.iupper(); ++i ) - { - EXPECT_EQ( x.get( i ), z.get( i ) ); - } - - // Testing scale, z = x - z.scale( 4.0 ); - for( globalIndex i = y.ilower(); i < y.iupper(); ++i ) - { - EXPECT_EQ( 4.0 * x.get( i ), z.get( i ) ); - } - - // Testing add/set single element - x.open(); - x.set( offset, -1 ); - x.close(); // set/add can't be interchanged - x.open(); - x.add( offset + 1, 10 ); - x.close(); - EXPECT_DOUBLE_EQ( x.get( offset ), -1 ); - EXPECT_DOUBLE_EQ( x.get( offset + 1 ), offset * 2 + 12 ); - - // Testing add/set c-style - { - globalIndex const inds[3] = { offset, offset + 1, offset + 2 }; - real64 const vals[3] = { -5.0, -6.0, 0.0 }; - y.zero(); - y.open(); - y.set( inds, vals, 2 ); - y.close(); - z.set( 1.0 ); - z.open(); - z.add( inds, vals, 2 ); - z.close(); - for( localIndex i = 0; i < 3; ++i ) - { - EXPECT_DOUBLE_EQ( y.get( inds[i] ), vals[i] ); - EXPECT_DOUBLE_EQ( z.get( inds[i] ), vals[i] + 1.0 ); - } - } - - // Testing add/set array1d-style - { - array1d< globalIndex > inds( 3 ); - inds[0] = offset; - inds[1] = offset + 1; - inds[2] = offset + 2; - array1d< real64 > vals( 3 ); - vals[0] = -5.0; - vals[1] = -6.0; - vals[2] = 0.0; - y.zero(); - y.open(); - y.set( inds, vals ); - y.close(); - z.set( 1.0 ); - z.open(); - z.add( inds, vals ); - z.close(); - for( localIndex i = 0; i < 3; ++i ) - { - EXPECT_DOUBLE_EQ( y.get( inds[i] ), vals[i] ); - EXPECT_DOUBLE_EQ( z.get( inds[i] ), vals[i] + 1.0 ); - } - } - - // Testing dot, axpy, axpby - x.set( 1.0 ); - y.set( 2.0 ); - z.set( 3.0 ); - - real64 const dotprod = x.dot( y ); - EXPECT_EQ( dotprod, 2 * y.globalSize() ); // sum_size 2 - - y.axpy( 2.0, x ); - for( globalIndex i = y.ilower(); i < y.iupper(); ++i ) - { - EXPECT_EQ( y.get( i ), 4.0 ); // 2*1 + 2 - } - - z.axpby( 2.0, x, 3.0 ); - for( globalIndex i = z.ilower(); i < z.iupper(); ++i ) - { - EXPECT_EQ( z.get( i ), 11.0 ); // 2*1 + 3*3 - } - - // Testing norms - x.zero(); - x.open(); - if( rank == 0 ) - { - globalIndex const inds2[2] = { 0, 1 }; - real64 const vals2[2] = { 3.0, -4.0 }; - x.set( inds2, vals2, 2 ); // 3, -4, 0 - } - x.close(); - EXPECT_EQ( x.norm1(), 7.0 ); - EXPECT_EQ( x.norm2(), 5.0 ); - EXPECT_EQ( x.normInf(), 4.0 ); - - // Testing extractLocalVector - real64 const * localVec = x.extractLocalVector(); - for( globalIndex i = x.ilower(); i < x.iupper(); ++i ) - { - EXPECT_EQ( localVec[x.getLocalRowID( i )], x.get( i ) ); - } -} - -#if 0 -TYPED_TEST_P( LAOperationsTest, MatrixFunctions ) -{ - // Define aliases - using Vector = typename TypeParam::ParallelVector; - using Matrix = typename TypeParam::ParallelMatrix; - - // Get the MPI rank - int numranks = MpiWrapper::Comm_size( MPI_COMM_GEOSX ); - int rank = MpiWrapper::Comm_rank( MPI_COMM_GEOSX ); - - std::cout << "*** Rank: " << rank << std::endl; - - // Dummy vector and Matrix - Matrix C; - Matrix D; - { - // Test matrix-matrix product: C = A*B - Matrix A; - compute2DLaplaceOperator( MPI_COMM_GEOSX, 2 * numranks, A ); - Matrix B( A ); - - A.multiply( B, C ); - A.leftMultiplyTranspose( A, D ); - } - - // Define some vectors, matrices - Vector vec1, vec2, vec3; - Matrix mat1, mat2, mat3, mat4; - mat1.createWithLocalSize( 2, 2, MPI_COMM_GEOSX ); // 2*numranks x 2*numranks - mat2.createWithGlobalSize( 2, 2, MPI_COMM_GEOSX ); // 2x2 - mat3.createWithLocalSize( 2, 3, 3, MPI_COMM_GEOSX ); // 2*numranks x 3*numranks - mat4.createWithGlobalSize( 3, 4, 3, MPI_COMM_GEOSX ); // 3x4 - - // Testing create, globalRows, globalCols - EXPECT_EQ( mat1.numGlobalRows(), 2 * numranks ); - EXPECT_EQ( mat1.numGlobalCols(), 2 * numranks ); - EXPECT_EQ( mat2.numGlobalRows(), 2 ); - EXPECT_EQ( mat2.numGlobalCols(), 2 ); - EXPECT_EQ( mat3.numGlobalRows(), 2 * numranks ); - EXPECT_EQ( mat3.numGlobalCols(), 3 * numranks ); - EXPECT_EQ( mat4.numGlobalRows(), 3 ); - EXPECT_EQ( mat4.numGlobalCols(), 4 ); - - // Testing add/set/insert element - // mat1.insert( 1, 0, .5 ); - // mat1.close(); - // mat1.set( 1, 0, 5 ); - // mat1.close(); - mat1.open(); - mat1.add( 1, 0, 1 ); - mat1.add( 1, 0, 2 ); - mat1.close(); - - // Testing add/set/insert c-style, getRowCopy - globalIndex inds1[2] = { 0, 2 }; - globalIndex inds2[1] = { 0 }; - globalIndex inds3[3] = { 0, 1, 2 }; - real64 vals1[2] = { 5, 10 }; - real64 vals2[1] = { 1 }; - real64 vals3[3] = { .5, 1, 2 }; - - globalIndex iRow = 1; - - if( ( mat4.ilower() <= iRow ) && ( iRow < mat4.iupper() ) ) - { - mat4.insert( iRow, inds3, vals3, 3 ); - // mat4.close(); - // mat4.open(); - mat4.set( iRow, inds1, vals1, 2 ); - // mat4.close(); - // mat4.open(); - mat4.add( iRow, inds2, vals2, 1 ); - // mat4.close(); - } - mat4.close(); - - array1d< real64 > colvals_CHECK( 3 ); - colvals_CHECK( 0 ) = 6; - colvals_CHECK( 1 ) = 1; - colvals_CHECK( 2 ) = 10; - - if( ( mat4.ilower() <= iRow ) && ( iRow < mat4.iupper() ) ) - { - localIndex const rowLength = mat.getGlobalRowLength( iRow ); - EXPECT_EQ( rowLength, colvals_CHECK.size() ); - array1d< real64 > colvals( rowLength ); - array1d< globalIndex > colinds( rowLength ); - mat4.getRowCopy( iRow, colinds, colvals ); - for( int i = 0; i < 3; ++i ) - { - EXPECT_DOUBLE_EQ( colvals( colinds[i] ), colvals_CHECK( i ) ); //HYPRE does not return sorted cols! - } - } - // Testing add/set/insert array1d - Matrix mat6; - mat6.createWithGlobalSize( 4, 4, MPI_COMM_GEOSX ); - array1d< real64 > vals6( 3 ); - array1d< real64 > vals7( 3 ); - array1d< globalIndex > inds6( 3 ); - vals6[0] = 1; - vals6[1] = .5; - vals6[2] = -3; - vals7[0] = 1; - vals7[1] = 1; - vals7[2] = 1; - inds6[0] = 0; - inds6[1] = 1; - inds6[2] = 3; - - iRow = 0; - if( ( mat6.ilower() <= iRow ) && ( iRow < mat6.iupper() ) ) - { - mat6.insert( iRow, inds6, vals6 ); - // mat6.close(); - mat6.set( iRow, inds6, vals7 ); - // mat6.close(); - mat6.add( iRow, inds6, vals6 ); - } - mat6.close(); - - // Testing add/set/insert array2d - Matrix mat7; - mat7.createWithGlobalSize( 4, 4, MPI_COMM_GEOSX ); - array1d< globalIndex > rows( 2 ); - array1d< globalIndex > cols( 2 ); - array2d< real64 > vals8( 2, 2 ); - rows[0] = 0; - rows[1] = 2; - cols[0] = 1; - cols[1] = 3; - vals8[0][0] = 1; - vals8[0][1] = 2; - vals8[1][0] = 3; - vals8[1][1] = 4; - if( ( mat7.ilower() <= *std::min_element( rows.data(), rows.data() + rows.size() ) ) && - ( *std::max_element( rows.data(), rows.data() + rows.size() ) < mat7.iupper() ) ) - { - mat7.insert( rows, cols, vals8 ); - // mat7.close(); - mat7.add( rows, cols, vals8 ); - // mat7.close(); - } - mat7.close(); - - // Testing set and zero - mat7.open(); - mat7.set( 2 ); - mat7.close(); - - mat7.open(); - mat1.zero(); - mat1.close(); - - // Testing vector multiply, matrix multiply, MatrixMatrixMultiply - vec1.createWithGlobalSize( 2, MPI_COMM_GEOSX ); - vec2.createWithGlobalSize( 2, MPI_COMM_GEOSX ); - vec1.set( 1 ); - vec1.close(); - globalIndex inds4[2] = { 0, 1 }; - real64 vals4[2] = { 1, 3 }; - real64 vals5[2] = { 2, 1 }; - - if( ( mat2.ilower() <= 0 ) && ( 0 < mat2.iupper() ) ) - { - mat2.insert( 0, inds4, vals4, 2 ); - } - if( ( mat2.ilower() <= 1 ) && ( 1 < mat2.iupper() ) ) - { - mat2.insert( 1, inds4, vals5, 2 ); - } - mat2.close(); - - mat2.multiply( vec1, vec2 ); - - if( ( vec2.ilower() <= 0 ) && ( 0 < vec2.iupper() ) ) - { - EXPECT_DOUBLE_EQ( vec2.get( 0 ), 4 ); - } - if( ( vec2.ilower() <= 1 ) && ( 1 < vec2.iupper() ) ) - { - EXPECT_DOUBLE_EQ( vec2.get( 1 ), 3 ); - } - - // Matrix-Matrix multiply - Matrix mat2mat2; - - { - Matrix mat22( mat2 ); - - mat22.multiply( mat22, mat2mat2 ); - } -} -#endif - -TYPED_TEST_P( LAOperationsTest, MatrixMatrixOperations ) -{ - using Matrix = typename TypeParam::ParallelMatrix; - - globalIndex const n = 100; - Matrix A; - compute2DLaplaceOperator( MPI_COMM_GEOSX, n, A ); - - Matrix A_squared; - A.multiply( A, A_squared ); - - real64 const a = A.normInf(); - real64 const b = A_squared.normInf(); - - EXPECT_DOUBLE_EQ( a, 8.0 ); - EXPECT_DOUBLE_EQ( b, 64.0 ); -} - -TYPED_TEST_P( LAOperationsTest, RectangularMatrixOperations ) -{ - using Matrix = typename TypeParam::ParallelMatrix; - - int mpiSize = MpiWrapper::Comm_size( MPI_COMM_GEOSX ); - - // Set a size that allows to run with arbitrary number of processes - globalIndex const nRows = std::max( 100, mpiSize ); - globalIndex const nCols = 2 * nRows; - - Matrix A; - A.createWithGlobalSize( nRows, nCols, 2, MPI_COMM_GEOSX ); - - A.open(); - for( globalIndex i = A.ilower(); i < A.iupper(); ++i ) - { - real64 const entry = static_cast< real64 >( i + 1 ); - A.insert( i, 2 * i, entry ); - A.insert( i, 2 * i + 1, -entry ); - } - A.close(); - - // Check on sizes - EXPECT_EQ( A.numGlobalRows(), nRows ); - EXPECT_EQ( A.numGlobalCols(), nCols ); - - // Check on norms - real64 const a = A.norm1(); - real64 const b = A.normInf(); - real64 const c = A.normFrobenius(); - - EXPECT_DOUBLE_EQ( a, static_cast< real64 >( nRows ) ); - EXPECT_DOUBLE_EQ( b, static_cast< real64 >( nCols ) ); - EXPECT_DOUBLE_EQ( c, std::sqrt( static_cast< real64 >( nRows * ( nRows + 1 ) * ( 2 * nRows + 1 ) ) / 3.0 ) ); -} - -REGISTER_TYPED_TEST_SUITE_P( LAOperationsTest, - VectorFunctions, - MatrixMatrixOperations, - RectangularMatrixOperations ); - -#ifdef GEOSX_USE_TRILINOS -INSTANTIATE_TYPED_TEST_SUITE_P( Trilinos, LAOperationsTest, TrilinosInterface, ); -#endif - -#ifdef GEOSX_USE_HYPRE -INSTANTIATE_TYPED_TEST_SUITE_P( Hypre, LAOperationsTest, HypreInterface, ); -#endif - -#ifdef GEOSX_USE_PETSC -INSTANTIATE_TYPED_TEST_SUITE_P( Petsc, LAOperationsTest, PetscInterface, ); -#endif - -/////////////////////////////////////////////////////////////////////////////////////// - -LinearSolverParameters params_Direct() -{ - LinearSolverParameters parameters; - parameters.solverType = geosx::LinearSolverParameters::SolverType::direct; - parameters.krylov.relTolerance = machinePrecision; - return parameters; -} - -LinearSolverParameters params_GMRES_ILU() -{ - LinearSolverParameters parameters; - parameters.krylov.relTolerance = 1e-8; - parameters.krylov.maxIterations = 300; - parameters.solverType = geosx::LinearSolverParameters::SolverType::gmres; - parameters.preconditionerType = geosx::LinearSolverParameters::PreconditionerType::iluk; - parameters.ilu.fill = 1; - return parameters; -} - -LinearSolverParameters params_GMRES_AMG() -{ - LinearSolverParameters parameters; - parameters.krylov.relTolerance = 1e-8; - parameters.krylov.maxIterations = 300; - parameters.solverType = geosx::LinearSolverParameters::SolverType::gmres; - parameters.preconditionerType = geosx::LinearSolverParameters::PreconditionerType::amg; - parameters.amg.smootherType = "gaussSeidel"; - parameters.amg.coarseType = "direct"; - return parameters; -} - -LinearSolverParameters params_CG_AMG() -{ - LinearSolverParameters parameters; - parameters.krylov.relTolerance = 1e-8; - parameters.krylov.maxIterations = 300; - parameters.solverType = geosx::LinearSolverParameters::SolverType::cg; - parameters.isSymmetric = true; - parameters.preconditionerType = geosx::LinearSolverParameters::PreconditionerType::amg; - parameters.amg.smootherType = "gaussSeidel"; - parameters.amg.coarseType = "direct"; - return parameters; -} - -template< typename LAI > -class SolverTestBase : public ::testing::Test -{ -public: - - using Matrix = typename LAI::ParallelMatrix; - using Vector = typename LAI::ParallelVector; - using Solver = typename LAI::LinearSolver; - -protected: - - Matrix matrix; - real64 cond_est = 1.0; - - void test( LinearSolverParameters const & params ) - { - // Create a random "true" solution vector - Vector sol_true; - sol_true.createWithGlobalSize( matrix.numGlobalCols(), matrix.getComm() ); - sol_true.rand(); - - // Create and compute the right-hand side vector - Vector rhs; - rhs.createWithGlobalSize( matrix.numGlobalRows(), matrix.getComm() ); - matrix.apply( sol_true, rhs ); - - // Create and zero out the computed solution vector - Vector sol_comp; - sol_comp.createWithGlobalSize( sol_true.globalSize(), sol_true.getComm() ); - sol_comp.zero(); - - // Create the solver and solve the system - Solver solver( params ); - solver.solve( matrix, sol_comp, rhs ); - EXPECT_TRUE( solver.result().success() ); - - // Check that solution is within epsilon of true - sol_comp.axpy( -1.0, sol_true ); - real64 const relTol = cond_est * params.krylov.relTolerance; - EXPECT_LT( sol_comp.norm2() / sol_true.norm2(), relTol ); - } -}; - -/////////////////////////////////////////////////////////////////////////////////////// - -template< typename LAI > -class SolverTestLaplace2D : public SolverTestBase< LAI > -{ -public: - - using Base = SolverTestBase< LAI >; - using Matrix = typename Base::Matrix; - using Vector = typename Base::Vector; - -protected: - - void SetUp() override - { - globalIndex constexpr n = 100; - compute2DLaplaceOperator( MPI_COMM_GEOSX, n, this->matrix ); - - // Condition number for the Laplacian matrix estimate: 4 * n^2 / pi^2 - this->cond_est = 4.0 * n * n / std::pow( M_PI, 2 ); - } -}; - -TYPED_TEST_SUITE_P( SolverTestLaplace2D ); - -TYPED_TEST_P( SolverTestLaplace2D, Direct ) -{ - this->test( params_Direct() ); -} - -TYPED_TEST_P( SolverTestLaplace2D, GMRES_ILU ) -{ - this->test( params_GMRES_ILU() ); -} - -TYPED_TEST_P( SolverTestLaplace2D, CG_AMG ) -{ - this->test( params_CG_AMG() ); -} - -REGISTER_TYPED_TEST_SUITE_P( SolverTestLaplace2D, - Direct, - GMRES_ILU, - CG_AMG ); - -#ifdef GEOSX_USE_TRILINOS -INSTANTIATE_TYPED_TEST_SUITE_P( Trilinos, SolverTestLaplace2D, TrilinosInterface, ); -#endif - -#ifdef GEOSX_USE_HYPRE -INSTANTIATE_TYPED_TEST_SUITE_P( Hypre, SolverTestLaplace2D, HypreInterface, ); -#endif - -#ifdef GEOSX_USE_PETSC -INSTANTIATE_TYPED_TEST_SUITE_P( Petsc, SolverTestLaplace2D, PetscInterface, ); -#endif - -/////////////////////////////////////////////////////////////////////////////////////// - -template< typename LAI > -class SolverTestElasticity2D : public SolverTestBase< LAI > -{ -public: - - using Base = SolverTestBase< LAI >; - using Matrix = typename Base::Matrix; - using Vector = typename Base::Vector; - -protected: - - void SetUp() override - { - globalIndex constexpr n = 100; - compute2DElasticityOperator( MPI_COMM_GEOSX, 1.0, 1.0, n, n, 10000., 0.2, this->matrix ); - - // Impose Dirichlet boundary conditions: fix domain bottom (first 2*(nCellsX + 1) rows of matrix) - this->matrix.open(); - for( globalIndex iRow = 0; iRow < 2 * (n + 1); ++iRow ) - { - if( this->matrix.getLocalRowID( iRow ) >= 0 ) - { - this->matrix.clearRow( iRow, true ); - } - } - this->matrix.close(); - this->cond_est = 1e3; // not a true condition number estimate, but enough to pass tests - } -}; - -TYPED_TEST_SUITE_P( SolverTestElasticity2D ); - -TYPED_TEST_P( SolverTestElasticity2D, Direct ) -{ - this->test( params_Direct() ); -} - -TYPED_TEST_P( SolverTestElasticity2D, GMRES_AMG ) -{ - LinearSolverParameters params = params_GMRES_AMG(); - params.amg.separateComponents = true; - params.dofsPerNode = 2; - this->test( params ); -} - -REGISTER_TYPED_TEST_SUITE_P( SolverTestElasticity2D, - Direct, - GMRES_AMG ); - -#ifdef GEOSX_USE_TRILINOS -INSTANTIATE_TYPED_TEST_SUITE_P( Trilinos, SolverTestElasticity2D, TrilinosInterface, ); -#endif - -#ifdef GEOSX_USE_HYPRE -INSTANTIATE_TYPED_TEST_SUITE_P( Hypre, SolverTestElasticity2D, HypreInterface, ); -#endif - -#ifdef GEOSX_USE_PETSC -INSTANTIATE_TYPED_TEST_SUITE_P( Petsc, SolverTestElasticity2D, PetscInterface, ); -#endif - -/////////////////////////////////////////////////////////////////////////////////////// -// -//template -//void parallel_vector_copy_constructor() -//{ -// -// int const rank = MpiWrapper::Comm_rank( MPI_COMM_WORLD ); -// typename LAI::ParallelVector x; -// typename LAI::ParallelVector y; -// localIndex const localSize = rank*10 + 1; -// -// array1d vec(localSize); -// -// // Populate vector with random coefficients -// BlasLapackLA::vectorRand( vec, -// BlasLapackLA::RandomNumberDistribution::UNIFORM_m1p1 ); -// -// x.create( vec, MPI_COMM_WORLD ); -// x.close(); -// -// y.create(x); -// -// for( globalIndex i = x.ilower(); i < x.iupper(); ++i ) -// { -// EXPECT_EQ( x.get( i ), y.get( i ) ); -// } -//} -// -//template -//void parallel_vector_create_with_local_size() -//{ -// int const rank = MpiWrapper::Comm_rank( MPI_COMM_WORLD ); -// typename LAI::ParallelVector x; -// localIndex const localSize = rank*10; -// globalIndex const globalSize = MpiWrapper::Sum( localSize ); -// -// x.createWithLocalSize( localSize, MPI_COMM_WORLD ); -// -// EXPECT_EQ( x.localSize(), localSize ); -// EXPECT_EQ( x.globalSize(), globalSize ); -//} -// -//template -//void parallel_vector_create_with_global_size() -//{ -// int const numRanks = MpiWrapper::Comm_size( MPI_COMM_WORLD ); -// typename LAI::ParallelVector x; -// localIndex const localSize = integer_conversion< localIndex >( numRanks ); -// globalIndex const globalSize = integer_conversion< globalIndex >( localSize * localSize ); -// -// // Testing createWithGlobalSize -// x.createWithGlobalSize( globalSize, MPI_COMM_WORLD ); -// EXPECT_EQ( x.localSize(), localSize ); -// EXPECT_EQ( x.globalSize(), globalSize ); -//} -// -//#ifdef GEOSX_USE_TRILINOS -// -//TEST( EpetraVector, EpetraVector ) -//{ -// parallel_vector_copy_constructor< TrilinosInterface >(); -//} -// -//TEST( EpetraVector, createWithLocalSize ) -//{ -// parallel_vector_create_with_local_size< TrilinosInterface >(); -//} -// -//TEST( EpetraVector, createWithGlobalSize ) -//{ -// parallel_vector_create_with_global_size< TrilinosInterface >(); -//} -// -//#endif -/////////////////////////////////////////////////////////////////////////////////////// - -int main( int argc, char * * argv ) -{ - ::testing::InitGoogleTest( &argc, argv ); - - geosx::basicSetup( argc, argv ); - - int const result = RUN_ALL_TESTS(); - geosx::basicCleanup(); - return result; -} diff --git a/src/coreComponents/linearAlgebra/unitTests/testLinearAlgebraUtils.hpp b/src/coreComponents/linearAlgebra/unitTests/testLinearAlgebraUtils.hpp index 723b87ccbd9..1685db19725 100644 --- a/src/coreComponents/linearAlgebra/unitTests/testLinearAlgebraUtils.hpp +++ b/src/coreComponents/linearAlgebra/unitTests/testLinearAlgebraUtils.hpp @@ -277,25 +277,29 @@ void compute2DElasticityOperator( MPI_Comm const comm, localIndex const rank = LvArray::integerConversion< localIndex >( MpiWrapper::Comm_rank( comm ) ); localIndex const nproc = LvArray::integerConversion< localIndex >( MpiWrapper::Comm_size( comm ) ); - // Compute total number of grid nodes (nNodes) and elements (nCells) - globalIndex const nCells = nCellsX * nCellsY; - GEOSX_ERROR_IF( nCells < nproc, "less than one cell per processor" ); - globalIndex const nNodes = ( nCellsX + 1 ) * ( nCellsY + 1); + GEOSX_ERROR_IF( nCellsY < nproc, "Less than one cell row per processor is not supported" ); + real64 const hx = domainSizeX / nCellsX; real64 const hy = domainSizeY / nCellsY; // Compute cell partitioning - localIndex const nLocalCells = LvArray::integerConversion< localIndex >( nCells / nproc ); + globalIndex const nCells = nCellsX * nCellsY; + localIndex const nLocalCells = LvArray::integerConversion< localIndex >( nCellsY / nproc * nCellsX ); localIndex const nExtraCells = LvArray::integerConversion< localIndex >( nCells ) - nLocalCells * nproc; globalIndex const iCellLower = rank * nLocalCells + ( rank == 0 ? 0 : nExtraCells ); - globalIndex const iCellUpper = iCellLower + nLocalCells + ( rank == 0 ? 0 : nExtraCells ) - 1; + globalIndex const iCellUpper = iCellLower + nLocalCells + ( rank == 0 ? 0 : nExtraCells ); + + // Compute node partitioning + globalIndex const iNodeLower = iCellLower / nCellsX + iCellLower; + globalIndex const iNodeUpper = iCellUpper / nCellsX + iCellUpper + ( rank == nproc - 1 ? nCellsX + 1 : 0 ); + localIndex const nLocalNodes = iNodeUpper - iNodeLower; // Construct local stiffness matrix (same for all cells) stackArray2d< real64, 8*8 > Ke( 8, 8 ); Q12d_local( hx, hy, youngModulus, poissonRatio, Ke ); // Create a matrix of global size N with at most 18 non-zeros per row - elasticity2D.createWithGlobalSize( nNodes*2, 18, comm ); + elasticity2D.createWithLocalSize( nLocalNodes * 2, 18, comm ); // Open the matrix elasticity2D.open(); @@ -304,11 +308,11 @@ void compute2DElasticityOperator( MPI_Comm const comm, stackArray1d< globalIndex, 4 > cellNodes( 4 ); stackArray1d< globalIndex, 8 > localDofIndex( 8 ); - for( localIndex iCell = iCellLower; iCell <= iCellUpper; ++iCell ) + for( globalIndex iCell = iCellLower; iCell < iCellUpper; ++iCell ) { // Compute local DOF global indeces - cellNodes( 0 ) = (iCell / nCellsX) + iCell; + cellNodes( 0 ) = iCell / nCellsX + iCell; cellNodes( 1 ) = cellNodes( 0 ) + 1; cellNodes( 3 ) = cellNodes( 1 ) + nCellsX; cellNodes( 2 ) = cellNodes( 3 ) + 1; diff --git a/src/coreComponents/linearAlgebra/unitTests/testMatrices.cpp b/src/coreComponents/linearAlgebra/unitTests/testMatrices.cpp new file mode 100644 index 00000000000..4b5dbb5663a --- /dev/null +++ b/src/coreComponents/linearAlgebra/unitTests/testMatrices.cpp @@ -0,0 +1,308 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2019 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2019 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2019 Total, S.A + * Copyright (c) 2019- GEOSX Contributors + * All right reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file testMatrices.cpp + */ + +#include + +#include "common/DataTypes.hpp" +#include "linearAlgebra/interfaces/InterfaceTypes.hpp" +#include "managers/initialization.hpp" + +#include "testLinearAlgebraUtils.hpp" + +using namespace geosx; + +template< typename LAI > +class LAOperationsTest : public ::testing::Test +{}; + +TYPED_TEST_SUITE_P( LAOperationsTest ); + +#if 0 +TYPED_TEST_P( LAOperationsTest, MatrixFunctions ) +{ + // Define aliases + using Vector = typename TypeParam::ParallelVector; + using Matrix = typename TypeParam::ParallelMatrix; + + // Get the MPI rank + int numranks = MpiWrapper::Comm_size( MPI_COMM_GEOSX ); + int rank = MpiWrapper::Comm_rank( MPI_COMM_GEOSX ); + + std::cout << "*** Rank: " << rank << std::endl; + + // Dummy vector and Matrix + Matrix C; + Matrix D; + { + // Test matrix-matrix product: C = A*B + Matrix A; + compute2DLaplaceOperator( MPI_COMM_GEOSX, 2 * numranks, A ); + Matrix B( A ); + + A.multiply( B, C ); + A.leftMultiplyTranspose( A, D ); + } + + // Define some vectors, matrices + Vector vec1, vec2, vec3; + Matrix mat1, mat2, mat3, mat4; + mat1.createWithLocalSize( 2, 2, MPI_COMM_GEOSX ); // 2*numranks x 2*numranks + mat2.createWithGlobalSize( 2, 2, MPI_COMM_GEOSX ); // 2x2 + mat3.createWithLocalSize( 2, 3, 3, MPI_COMM_GEOSX ); // 2*numranks x 3*numranks + mat4.createWithGlobalSize( 3, 4, 3, MPI_COMM_GEOSX ); // 3x4 + + // Testing create, globalRows, globalCols + EXPECT_EQ( mat1.numGlobalRows(), 2 * numranks ); + EXPECT_EQ( mat1.numGlobalCols(), 2 * numranks ); + EXPECT_EQ( mat2.numGlobalRows(), 2 ); + EXPECT_EQ( mat2.numGlobalCols(), 2 ); + EXPECT_EQ( mat3.numGlobalRows(), 2 * numranks ); + EXPECT_EQ( mat3.numGlobalCols(), 3 * numranks ); + EXPECT_EQ( mat4.numGlobalRows(), 3 ); + EXPECT_EQ( mat4.numGlobalCols(), 4 ); + + // Testing add/set/insert element + // mat1.insert( 1, 0, .5 ); + // mat1.close(); + // mat1.set( 1, 0, 5 ); + // mat1.close(); + mat1.open(); + mat1.add( 1, 0, 1 ); + mat1.add( 1, 0, 2 ); + mat1.close(); + + // Testing add/set/insert c-style, getRowCopy + globalIndex inds1[2] = { 0, 2 }; + globalIndex inds2[1] = { 0 }; + globalIndex inds3[3] = { 0, 1, 2 }; + real64 vals1[2] = { 5, 10 }; + real64 vals2[1] = { 1 }; + real64 vals3[3] = { .5, 1, 2 }; + + globalIndex iRow = 1; + + if( ( mat4.ilower() <= iRow ) && ( iRow < mat4.iupper() ) ) + { + mat4.insert( iRow, inds3, vals3, 3 ); + // mat4.close(); + // mat4.open(); + mat4.set( iRow, inds1, vals1, 2 ); + // mat4.close(); + // mat4.open(); + mat4.add( iRow, inds2, vals2, 1 ); + // mat4.close(); + } + mat4.close(); + + array1d< real64 > colvals_CHECK( 3 ); + colvals_CHECK( 0 ) = 6; + colvals_CHECK( 1 ) = 1; + colvals_CHECK( 2 ) = 10; + + if( ( mat4.ilower() <= iRow ) && ( iRow < mat4.iupper() ) ) + { + localIndex const rowLength = mat.getGlobalRowLength( iRow ); + EXPECT_EQ( rowLength, colvals_CHECK.size() ); + array1d< real64 > colvals( rowLength ); + array1d< globalIndex > colinds( rowLength ); + mat4.getRowCopy( iRow, colinds, colvals ); + for( int i = 0; i < 3; ++i ) + { + EXPECT_DOUBLE_EQ( colvals( colinds[i] ), colvals_CHECK( i ) ); //HYPRE does not return sorted cols! + } + } + // Testing add/set/insert array1d + Matrix mat6; + mat6.createWithGlobalSize( 4, 4, MPI_COMM_GEOSX ); + array1d< real64 > vals6( 3 ); + array1d< real64 > vals7( 3 ); + array1d< globalIndex > inds6( 3 ); + vals6[0] = 1; + vals6[1] = .5; + vals6[2] = -3; + vals7[0] = 1; + vals7[1] = 1; + vals7[2] = 1; + inds6[0] = 0; + inds6[1] = 1; + inds6[2] = 3; + + iRow = 0; + if( ( mat6.ilower() <= iRow ) && ( iRow < mat6.iupper() ) ) + { + mat6.insert( iRow, inds6, vals6 ); + // mat6.close(); + mat6.set( iRow, inds6, vals7 ); + // mat6.close(); + mat6.add( iRow, inds6, vals6 ); + } + mat6.close(); + + // Testing add/set/insert array2d + Matrix mat7; + mat7.createWithGlobalSize( 4, 4, MPI_COMM_GEOSX ); + array1d< globalIndex > rows( 2 ); + array1d< globalIndex > cols( 2 ); + array2d< real64 > vals8( 2, 2 ); + rows[0] = 0; + rows[1] = 2; + cols[0] = 1; + cols[1] = 3; + vals8[0][0] = 1; + vals8[0][1] = 2; + vals8[1][0] = 3; + vals8[1][1] = 4; + if( ( mat7.ilower() <= *std::min_element( rows.data(), rows.data() + rows.size() ) ) && + ( *std::max_element( rows.data(), rows.data() + rows.size() ) < mat7.iupper() ) ) + { + mat7.insert( rows, cols, vals8 ); + // mat7.close(); + mat7.add( rows, cols, vals8 ); + // mat7.close(); + } + mat7.close(); + + // Testing set and zero + mat7.open(); + mat7.set( 2 ); + mat7.close(); + + mat7.open(); + mat1.zero(); + mat1.close(); + + // Testing vector multiply, matrix multiply, MatrixMatrixMultiply + vec1.createWithGlobalSize( 2, MPI_COMM_GEOSX ); + vec2.createWithGlobalSize( 2, MPI_COMM_GEOSX ); + vec1.set( 1 ); + vec1.close(); + globalIndex inds4[2] = { 0, 1 }; + real64 vals4[2] = { 1, 3 }; + real64 vals5[2] = { 2, 1 }; + + if( ( mat2.ilower() <= 0 ) && ( 0 < mat2.iupper() ) ) + { + mat2.insert( 0, inds4, vals4, 2 ); + } + if( ( mat2.ilower() <= 1 ) && ( 1 < mat2.iupper() ) ) + { + mat2.insert( 1, inds4, vals5, 2 ); + } + mat2.close(); + + mat2.multiply( vec1, vec2 ); + + if( ( vec2.ilower() <= 0 ) && ( 0 < vec2.iupper() ) ) + { + EXPECT_DOUBLE_EQ( vec2.get( 0 ), 4 ); + } + if( ( vec2.ilower() <= 1 ) && ( 1 < vec2.iupper() ) ) + { + EXPECT_DOUBLE_EQ( vec2.get( 1 ), 3 ); + } + + // Matrix-Matrix multiply + Matrix mat2mat2; + + { + Matrix mat22( mat2 ); + + mat22.multiply( mat22, mat2mat2 ); + } +} +#endif + +TYPED_TEST_P( LAOperationsTest, MatrixMatrixOperations ) +{ + using Matrix = typename TypeParam::ParallelMatrix; + + globalIndex const n = 100; + Matrix A; + compute2DLaplaceOperator( MPI_COMM_GEOSX, n, A ); + + Matrix A_squared; + A.multiply( A, A_squared ); + + real64 const a = A.normInf(); + real64 const b = A_squared.normInf(); + + EXPECT_DOUBLE_EQ( a, 8.0 ); + EXPECT_DOUBLE_EQ( b, 64.0 ); +} + +TYPED_TEST_P( LAOperationsTest, RectangularMatrixOperations ) +{ + using Matrix = typename TypeParam::ParallelMatrix; + + int mpiSize = MpiWrapper::Comm_size( MPI_COMM_GEOSX ); + + // Set a size that allows to run with arbitrary number of processes + globalIndex const nRows = std::max( 100, mpiSize ); + globalIndex const nCols = 2 * nRows; + + Matrix A; + A.createWithGlobalSize( nRows, nCols, 2, MPI_COMM_GEOSX ); + + A.open(); + for( globalIndex i = A.ilower(); i < A.iupper(); ++i ) + { + real64 const entry = static_cast< real64 >( i + 1 ); + A.insert( i, 2 * i, entry ); + A.insert( i, 2 * i + 1, -entry ); + } + A.close(); + + // Check on sizes + EXPECT_EQ( A.numGlobalRows(), nRows ); + EXPECT_EQ( A.numGlobalCols(), nCols ); + + // Check on norms + real64 const a = A.norm1(); + real64 const b = A.normInf(); + real64 const c = A.normFrobenius(); + + EXPECT_DOUBLE_EQ( a, static_cast< real64 >( nRows ) ); + EXPECT_DOUBLE_EQ( b, static_cast< real64 >( nCols ) ); + EXPECT_DOUBLE_EQ( c, std::sqrt( static_cast< real64 >( nRows * ( nRows + 1 ) * ( 2 * nRows + 1 ) ) / 3.0 ) ); +} + +REGISTER_TYPED_TEST_SUITE_P( LAOperationsTest, + MatrixMatrixOperations, + RectangularMatrixOperations ); + +#ifdef GEOSX_USE_TRILINOS +INSTANTIATE_TYPED_TEST_SUITE_P( Trilinos, LAOperationsTest, TrilinosInterface, ); +INSTANTIATE_TYPED_TEST_SUITE_P( TrilinosTpetra, LAOperationsTest, TrilinosTpetraInterface, ); +#endif + +#ifdef GEOSX_USE_HYPRE +INSTANTIATE_TYPED_TEST_SUITE_P( Hypre, LAOperationsTest, HypreInterface, ); +#endif + +#ifdef GEOSX_USE_PETSC +INSTANTIATE_TYPED_TEST_SUITE_P( Petsc, LAOperationsTest, PetscInterface, ); +#endif + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + geosx::basicSetup( argc, argv ); + int const result = RUN_ALL_TESTS(); + geosx::basicCleanup(); + return result; +} diff --git a/src/coreComponents/linearAlgebra/unitTests/testVectors.cpp b/src/coreComponents/linearAlgebra/unitTests/testVectors.cpp new file mode 100644 index 00000000000..6624428edce --- /dev/null +++ b/src/coreComponents/linearAlgebra/unitTests/testVectors.cpp @@ -0,0 +1,601 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2019 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2019 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2019 Total, S.A + * Copyright (c) 2019- GEOSX Contributors + * All right reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file testVectors.cpp + */ + +#include + +#include "common/DataTypes.hpp" +#include "linearAlgebra/interfaces/InterfaceTypes.hpp" +#include "managers/initialization.hpp" + +using namespace geosx; + +/** ---------------------- Helpers ---------------------- **/ + +array1d< real64 > makeLocalValues( localIndex const localSize, globalIndex const rankOffset ) +{ + array1d< real64 > localValues( localSize ); + for( localIndex i = 0; i < localSize; ++i ) + { + globalIndex const row = rankOffset + i; + localValues[i] = std::pow( -1.0, row ) * ( 1.0 + row ); + } + return localValues; +} + +array1d< real64 > makeLocalValuesUniform( localIndex const localSize ) +{ + int const rank = MpiWrapper::Comm_rank( MPI_COMM_GEOSX ); + localIndex const rankOffset = localSize * rank; + + return makeLocalValues( localSize, rankOffset ); +} + +array1d< real64 > makeLocalValuesNonUniform( localIndex const startSize ) +{ + int const rank = MpiWrapper::Comm_rank( MPI_COMM_GEOSX ); + localIndex const localSize = rank + startSize; + globalIndex const rankOffset = ( rank + startSize ) * ( rank + startSize - 1 ) / 2; + + return makeLocalValues( localSize, rankOffset ); +} + +namespace ops +{ + +auto identity = []( real64 const x ){ return x; }; + +struct multiply +{ + real64 factor; + real64 operator()( real64 const x ) const { return factor * x; } +}; + +auto reciprocal = []( real64 const x ){ return 1.0 / x; }; + +} + +template< typename OP_X = decltype( ops::identity ), typename OP_Y = decltype( ops::identity ) > +void compareValues( arrayView1d< real64 > const & x, + arrayView1d< real64 > const & y, + bool exact = true, + OP_X const op_x = ops::identity, + OP_Y const op_y = ops::identity ) +{ + EXPECT_EQ( x.size(), y.size() ); + x.move( LvArray::MemorySpace::CPU, false ); + y.move( LvArray::MemorySpace::CPU, false ); + for( localIndex i = 0; i < x.size(); ++i ) + { + if( exact ) + { + EXPECT_EQ( op_x( x[i] ), op_y( y[i] ) ); + } + else + { + EXPECT_DOUBLE_EQ( op_x( x[i] ), op_y( y[i] ) ); + } + } +} + +template< typename VEC, typename OP_X = decltype( ops::identity ), typename OP_Y = decltype( ops::identity ) > +void compareVectors( VEC const & x, + VEC const & y, + bool exact = true, + OP_X const op_x = ops::identity, + OP_Y const op_y = ops::identity ) +{ + EXPECT_TRUE( MpiWrapper::Comm_compare( y.getComm(), x.getComm() ) ); + EXPECT_EQ( y.localSize(), x.localSize() ); + EXPECT_EQ( y.globalSize(), x.globalSize() ); + + array1d< real64 > xval( x.localSize() ); + x.extract( xval ); + + array1d< real64 > yval( y.localSize() ); + y.extract( yval ); + + compareValues( xval, yval, exact, op_x, op_y ); +} + +/** ---------------------- Tests ---------------------- **/ + +template< typename LAI > +class VectorTest : public ::testing::Test +{}; + +TYPED_TEST_SUITE_P( VectorTest ); + +TYPED_TEST_P( VectorTest, createWithLocalSize ) +{ + using Vector = typename TypeParam::ParallelVector; + + MPI_Comm const comm = MPI_COMM_GEOSX; + int const rank = MpiWrapper::Comm_rank( comm ); + int const nproc = MpiWrapper::Comm_size( comm ); + + localIndex const localSize = rank + 1; + globalIndex const globalSize = ( nproc + 1 ) * nproc / 2; + globalIndex const rankOffset = ( rank + 1 ) * rank / 2; + + Vector x; + x.createWithLocalSize( localSize, comm ); + + EXPECT_TRUE( x.ready() ); + EXPECT_TRUE( MpiWrapper::Comm_compare( x.getComm(), comm ) ); + EXPECT_EQ( x.localSize(), localSize ); + EXPECT_EQ( x.globalSize(), globalSize ); + EXPECT_EQ( x.ilower(), rankOffset ); + EXPECT_EQ( x.iupper(), rankOffset + localSize ); +} + +TYPED_TEST_P( VectorTest, createWithGlobalSize ) +{ + using Vector = typename TypeParam::ParallelVector; + + MPI_Comm const comm = MPI_COMM_GEOSX; + int const rank = MpiWrapper::Comm_rank( comm ); + int const nproc = MpiWrapper::Comm_size( comm ); + + localIndex const localSize = 3; + globalIndex const globalSize = localSize * nproc; + globalIndex const rankOffset = localSize * rank; + + Vector x; + x.createWithGlobalSize( globalSize, comm ); + + EXPECT_TRUE( x.ready() ); + EXPECT_TRUE( MpiWrapper::Comm_compare( x.getComm(), comm ) ); + EXPECT_EQ( x.localSize(), localSize ); + EXPECT_EQ( x.globalSize(), globalSize ); + EXPECT_EQ( x.ilower(), rankOffset ); + EXPECT_EQ( x.iupper(), rankOffset + localSize ); +} + +TYPED_TEST_P( VectorTest, createWithLocalValues ) +{ + using Vector = typename TypeParam::ParallelVector; + + array1d< real64 > const valuesInitial = makeLocalValuesNonUniform( 1 ); + localIndex const localSize = valuesInitial.size(); + globalIndex const globalSize = MpiWrapper::Sum( localSize, MPI_COMM_GEOSX ); + + Vector x; + x.createWithLocalValues( valuesInitial, MPI_COMM_GEOSX ); + + EXPECT_TRUE( x.ready() ); + EXPECT_TRUE( MpiWrapper::Comm_compare( x.getComm(), MPI_COMM_GEOSX ) ); + EXPECT_EQ( x.localSize(), localSize ); + EXPECT_EQ( x.globalSize(), globalSize ); + + array1d< real64 > const valuesExtracted( localSize ); + x.extract( valuesExtracted ); + + compareValues( valuesExtracted, valuesInitial ); +} + +TYPED_TEST_P( VectorTest, copyConstruction ) +{ + using Vector = typename TypeParam::ParallelVector; + + array1d< real64 > const localValues = makeLocalValuesUniform( 3 ); + + Vector x; + x.createWithLocalValues( localValues, MPI_COMM_GEOSX ); + + Vector y( x ); + + // Test that values are equal after copy contruction + EXPECT_TRUE( y.ready() ); + compareVectors( x, y ); + + // Test that vectors don't share data after copy construction + real64 const factor = 0.5; + x.scale( factor ); + + compareVectors( x, y, true, ops::identity, ops::multiply{factor} ); +} + +TYPED_TEST_P( VectorTest, moveConstruction ) +{ + using Vector = typename TypeParam::ParallelVector; + + array1d< real64 > const valuesInitial = makeLocalValuesUniform( 3 ); + localIndex const localSize = valuesInitial.size(); + globalIndex const globalSize = MpiWrapper::Sum( localSize ); + + Vector x; + x.createWithLocalValues( valuesInitial, MPI_COMM_GEOSX ); + + Vector y( std::move( x ) ); + + EXPECT_TRUE( y.ready() ); + EXPECT_TRUE( MpiWrapper::Comm_compare( y.getComm(), MPI_COMM_GEOSX ) ); + EXPECT_EQ( y.localSize(), localSize ); + EXPECT_EQ( y.globalSize(), globalSize ); + + array1d< real64 > const valuesExtracted( localSize ); + y.extract( valuesExtracted ); + + compareValues( valuesExtracted, valuesInitial ); +} + +TYPED_TEST_P( VectorTest, copy ) +{ + using Vector = typename TypeParam::ParallelVector; + + localIndex const localSize = 3; + + Vector x; + x.createWithLocalSize( localSize, MPI_COMM_GEOSX ); + x.rand(); + + Vector y; + y.createWithLocalSize( localSize, MPI_COMM_GEOSX ); + y.copy( x ); + + compareVectors( x, y ); +} + +TYPED_TEST_P( VectorTest, extract ) +{ + using Vector = typename TypeParam::ParallelVector; + + array1d< real64 > const valuesInitial = makeLocalValuesUniform( 3 ); + + Vector x; + x.createWithLocalValues( valuesInitial, MPI_COMM_GEOSX ); + + array1d< real64 > const valuesExtracted( valuesInitial.size() ); + x.extract( valuesExtracted ); + + compareValues( valuesExtracted, valuesInitial ); +} + +TYPED_TEST_P( VectorTest, localGlobalRowID ) +{ + using Vector = typename TypeParam::ParallelVector; + + int const rank = MpiWrapper::Comm_rank( MPI_COMM_GEOSX ); + + localIndex const localSize = rank + 1; + globalIndex const offset = ( rank + 1 ) * rank / 2; + + Vector x; + x.createWithLocalSize( localSize, MPI_COMM_GEOSX ); + + for( localIndex i = 0; i < localSize; ++i ) + { + EXPECT_EQ( x.getGlobalRowID( i ), offset + i ); + EXPECT_EQ( x.getLocalRowID( offset + i ), i ); + } +} + +TYPED_TEST_P( VectorTest, getSingleValue ) +{ + using Vector = typename TypeParam::ParallelVector; + + localIndex const localSize = 3; + array1d< real64 > const localValues = makeLocalValuesUniform( 3 ); + array1d< real64 > const localValuesCopy = localValues; + + Vector x; + x.createWithLocalValues( localValues, MPI_COMM_GEOSX ); + + for( localIndex i = 0; i < localSize; ++i ) + { + EXPECT_EQ( x.get( x.getGlobalRowID( i ) ), localValuesCopy[i] ); + } +} + +TYPED_TEST_P( VectorTest, getMultipleValues ) +{ + using Vector = typename TypeParam::ParallelVector; + + localIndex const localSize = 3; + array1d< real64 > const valuesInitial = makeLocalValuesUniform( 3 ); + + Vector x; + x.createWithLocalValues( valuesInitial, MPI_COMM_GEOSX ); + + array1d< globalIndex > const rowIndices( localSize ); + array1d< real64 > const valuesExtracted( localSize ); + + for( localIndex i = 0; i < localSize; ++i ) + { + rowIndices[i] = x.getGlobalRowID( i ); + } + + x.get( rowIndices, valuesExtracted ); + compareValues( valuesExtracted, valuesInitial ); +} + +TYPED_TEST_P( VectorTest, setSingleValue ) +{ + using Vector = typename TypeParam::ParallelVector; + + localIndex const localSize = 3; + + Vector x; + x.createWithLocalSize( localSize, MPI_COMM_GEOSX ); + + x.open(); + for( localIndex i = 0; i < localSize; ++i ) + { + x.set( x.getGlobalRowID( i ), x.ilower() + i ); + } + x.close(); + + for( localIndex i = 0; i < localSize; ++i ) + { + EXPECT_EQ( x.get( x.getGlobalRowID( i ) ), x.ilower() + i ); + } +} + +TYPED_TEST_P( VectorTest, addSingleValue ) +{ + using Vector = typename TypeParam::ParallelVector; + + localIndex const localSize = 3; + + Vector x; + x.createWithLocalSize( localSize, MPI_COMM_GEOSX ); + x.set( 1.0 ); + + x.open(); + for( localIndex i = 0; i < localSize; ++i ) + { + x.add( x.getGlobalRowID( i ), x.ilower() + i ); + } + x.close(); + + for( localIndex i = 0; i < localSize; ++i ) + { + EXPECT_EQ( x.get( x.getGlobalRowID( i ) ), 1.0 + x.ilower() + i ); + } +} + +TYPED_TEST_P( VectorTest, setAllValues ) +{ + using Vector = typename TypeParam::ParallelVector; + + localIndex const localSize = 3; + real64 const value = 1.23; + + Vector x; + x.createWithLocalSize( localSize, MPI_COMM_GEOSX ); + x.set( value ); + + array1d< real64 > const valuesExtracted( localSize ); + x.extract( valuesExtracted ); + valuesExtracted.move( LvArray::MemorySpace::CPU, false ); + + for( localIndex i = 0; i < localSize; ++i ) + { + EXPECT_EQ( valuesExtracted[i], value ); + } +} + +TYPED_TEST_P( VectorTest, zeroAllValues ) +{ + using Vector = typename TypeParam::ParallelVector; + + localIndex const localSize = 3; + + Vector x; + x.createWithLocalSize( localSize, MPI_COMM_GEOSX ); + x.zero(); + + array1d< real64 > const valuesExtracted( localSize ); + x.extract( valuesExtracted ); + valuesExtracted.move( LvArray::MemorySpace::CPU, false ); + + for( localIndex i = 0; i < localSize; ++i ) + { + EXPECT_EQ( valuesExtracted[i], 0.0 ); + } +} + +TYPED_TEST_P( VectorTest, scaleValues ) +{ + using Vector = typename TypeParam::ParallelVector; + + localIndex const localSize = 3; + array1d< real64 > const valuesInitial = makeLocalValuesUniform( 3 ); + array1d< real64 > const valuesCopy = valuesInitial; + + Vector x; + x.createWithLocalValues( valuesInitial, MPI_COMM_GEOSX ); + + real64 const factor = 0.5; + x.scale( factor ); + + array1d< real64 > const valuesExtracted( localSize ); + x.extract( valuesExtracted ); + + compareValues( valuesExtracted, valuesCopy, true, ops::identity, ops::multiply{factor} ); +} + +TYPED_TEST_P( VectorTest, reciprocal ) +{ + using Vector = typename TypeParam::ParallelVector; + + localIndex const localSize = 3; + array1d< real64 > const valuesInitial = makeLocalValuesUniform( 3 ); + array1d< real64 > const valuesCopy = valuesInitial; + + Vector x; + x.createWithLocalValues( valuesInitial, MPI_COMM_GEOSX ); + + x.reciprocal(); + + array1d< real64 > const valuesExtracted( localSize ); + x.extract( valuesExtracted ); + + compareValues( valuesExtracted, valuesCopy, true, ops::reciprocal ); +} + +TYPED_TEST_P( VectorTest, dotProduct ) +{ + using Vector = typename TypeParam::ParallelVector; + + array1d< real64 > const xValues = makeLocalValuesUniform( 3 ); + array1d< real64 > const yValues = xValues; + + Vector x; + x.createWithLocalValues( xValues, MPI_COMM_GEOSX ); + + Vector y; + y.createWithLocalValues( yValues, MPI_COMM_GEOSX ); + y.reciprocal(); + + real64 const dp = x.dot( y ); + EXPECT_DOUBLE_EQ( dp, x.globalSize() ); +} + +TYPED_TEST_P( VectorTest, axpy ) +{ + using Vector = typename TypeParam::ParallelVector; + + real64 const alpha = 2.0; + + array1d< real64 > const xValues = makeLocalValuesUniform( 3 ); + array1d< real64 > const yValues = xValues; + + Vector x; + x.createWithLocalValues( xValues, MPI_COMM_GEOSX ); + + Vector y; + y.createWithLocalValues( yValues, MPI_COMM_GEOSX ); + + x.axpy( alpha, y ); + + compareVectors( y, x, false, ops::multiply{ 1.0 + alpha } ); +} + +TYPED_TEST_P( VectorTest, axpby ) +{ + using Vector = typename TypeParam::ParallelVector; + + real64 const alpha = 2.0; + real64 const beta = 3.0; + + array1d< real64 > const xValues = makeLocalValuesUniform( 3 ); + array1d< real64 > const yValues = xValues; + + Vector x; + x.createWithLocalValues( xValues, MPI_COMM_GEOSX ); + + Vector y; + y.createWithLocalValues( yValues, MPI_COMM_GEOSX ); + + x.axpby( alpha, y, beta ); + + compareVectors( y, x, false, ops::multiply{ alpha + beta } ); +} + +TYPED_TEST_P( VectorTest, norm1 ) +{ + using Vector = typename TypeParam::ParallelVector; + + array1d< real64 > const xValues = makeLocalValuesUniform( 3 ); + + Vector x; + x.createWithLocalValues( xValues, MPI_COMM_GEOSX ); + x.scale( -1.0 ); + + real64 const normTrue = ( x.globalSize() + 1 ) * x.globalSize() / 2; + + EXPECT_DOUBLE_EQ( x.norm1(), normTrue ); +} + +TYPED_TEST_P( VectorTest, norm2 ) +{ + using Vector = typename TypeParam::ParallelVector; + + array1d< real64 > const xValues = makeLocalValuesUniform( 3 ); + + Vector x; + x.createWithLocalValues( xValues, MPI_COMM_GEOSX ); + x.scale( -1.0 ); + + real64 const normTrue = std::sqrt( ( 2 * x.globalSize() + 1 ) * ( x.globalSize() + 1 ) * x.globalSize() / 6 ); + + EXPECT_DOUBLE_EQ( x.norm2(), normTrue ); +} + +TYPED_TEST_P( VectorTest, normInf ) +{ + using Vector = typename TypeParam::ParallelVector; + + array1d< real64 > const xValues = makeLocalValuesUniform( 3 ); + + Vector x; + x.createWithLocalValues( xValues, MPI_COMM_GEOSX ); + x.scale( -1.0 ); + + real64 const normTrue = x.globalSize(); + + EXPECT_DOUBLE_EQ( x.normInf(), normTrue ); +} + +REGISTER_TYPED_TEST_SUITE_P( VectorTest, + createWithLocalSize, + createWithGlobalSize, + createWithLocalValues, + copyConstruction, + moveConstruction, + copy, + extract, + localGlobalRowID, + getSingleValue, + getMultipleValues, + setSingleValue, + addSingleValue, + setAllValues, + zeroAllValues, + scaleValues, + reciprocal, + dotProduct, + axpy, + axpby, + norm1, + norm2, + normInf ); + +#ifdef GEOSX_USE_TRILINOS +INSTANTIATE_TYPED_TEST_SUITE_P( Trilinos, VectorTest, TrilinosInterface, ); +INSTANTIATE_TYPED_TEST_SUITE_P( TrilinosTpetra, VectorTest, TrilinosTpetraInterface, ); +#endif + +#ifdef GEOSX_USE_HYPRE +INSTANTIATE_TYPED_TEST_SUITE_P( Hypre, VectorTest, HypreInterface, ); +#endif + +#ifdef GEOSX_USE_PETSC +INSTANTIATE_TYPED_TEST_SUITE_P( Petsc, VectorTest, PetscInterface, ); +#endif + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + geosx::basicSetup( argc, argv ); + chai::ArrayManager::getInstance()->disableCallbacks(); + int const result = RUN_ALL_TESTS(); + geosx::basicCleanup(); + return result; +} diff --git a/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp b/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp index ba31e22fb2d..2cf607cf44f 100644 --- a/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp +++ b/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp @@ -50,17 +50,19 @@ struct LinearSolverParameters */ enum class PreconditionerType : integer { - none, ///< No preconditioner - jacobi, ///< Jacobi smoothing - gs, ///< Gauss-Seidel smoothing - sgs, ///< Symmetric Gauss-Seidel smoothing - iluk, ///< Incomplete LU with k-level of fill - ilut, ///< Incomplete LU with thresholding - icc, ///< Incomplete Cholesky - ict, ///< Incomplete Cholesky with thresholding - amg, ///< Algebraic Multigrid - mgr, ///< Multigrid reduction (Hypre only) - block ///< Block preconditioner + none, ///< No preconditioner + jacobi, ///< Jacobi smoothing + gs, ///< Gauss-Seidel smoothing + sgs, ///< Symmetric Gauss-Seidel smoothing + chebyshev, ///< Chebyshev polynomial smoothing + iluk, ///< Incomplete LU with k-level of fill + ilut, ///< Incomplete LU with thresholding + icc, ///< Incomplete Cholesky + ict, ///< Incomplete Cholesky with thresholding + amg, ///< Algebraic Multigrid + mgr, ///< Multigrid reduction (Hypre only) + block, ///< Block preconditioner + direct, ///< Direct solver as preconditioner }; integer logLevel = 0; ///< Output level [0=none, 1=basic, 2=everything] @@ -128,14 +130,14 @@ struct LinearSolverParameters /// Algebraic multigrid parameters struct AMG { + PreconditionerType smootherType = PreconditionerType::gs; ///< Smoother type + PreconditionerType coarseType = PreconditionerType::direct; ///< Coarse-level solver/smoother + integer maxLevels = 20; ///< Maximum number of coarsening levels string cycleType = "V"; ///< AMG cycle type - string smootherType = "gaussSeidel"; ///< Smoother type - string coarseType = "direct"; ///< Coarse-level solver/smoother integer numSweeps = 2; ///< Number of smoother sweeps string preOrPostSmoothing = "both"; ///< Pre and/or post smoothing [pre,post,both] - real64 threshold = 0.0; ///< Threshold for "strong connections" (for classical and - ///< smoothed-aggregation AMG) + real64 threshold = 0.0; ///< Threshold for "strong connections" (for classical and smoothed-aggregation AMG) integer separateComponents = false; ///< Apply a separate component filter before AMG construction string nullSpaceType = "constantModes"; ///< Null space type [constantModes,rigidBodyModes] } @@ -179,13 +181,15 @@ ENUM_STRINGS( LinearSolverParameters::PreconditionerType, "jacobi", "gs", "sgs", + "chebyshev", "iluk", "ilut", "icc", "ict", "amg", "mgr", - "block" ) + "block", + "direct" ) ENUM_STRINGS( LinearSolverParameters::Direct::ColPerm, "none", diff --git a/src/coreComponents/mpiCommunications/MpiWrapper.hpp b/src/coreComponents/mpiCommunications/MpiWrapper.hpp index 1ef95c8ca7f..a9046682579 100644 --- a/src/coreComponents/mpiCommunications/MpiWrapper.hpp +++ b/src/coreComponents/mpiCommunications/MpiWrapper.hpp @@ -178,6 +178,17 @@ struct MpiWrapper return size; } + inline static bool Comm_compare( MPI_Comm const & comm1, MPI_Comm const & comm2 ) + { +#ifdef GEOSX_USE_MPI + int result; + MPI_Comm_compare( comm1, comm2, &result ); + return result == MPI_IDENT || result == MPI_CONGRUENT; +#else + return comm1 == comm2; +#endif + } + static int Init( int * argc, char * * * argv ); static void Finalize(); diff --git a/src/coreComponents/physicsSolvers/GEOSX_PTP b/src/coreComponents/physicsSolvers/GEOSX_PTP index d1f2a7cbeab..f67d3e5212c 160000 --- a/src/coreComponents/physicsSolvers/GEOSX_PTP +++ b/src/coreComponents/physicsSolvers/GEOSX_PTP @@ -1 +1 @@ -Subproject commit d1f2a7cbeabd7c205521654942429c3255d99483 +Subproject commit f67d3e5212cb06e59e555008faae8669572c49f2 diff --git a/src/coreComponents/physicsSolvers/LinearSolverParameters.cpp b/src/coreComponents/physicsSolvers/LinearSolverParameters.cpp index 36a66097f4d..f0cd664913e 100644 --- a/src/coreComponents/physicsSolvers/LinearSolverParameters.cpp +++ b/src/coreComponents/physicsSolvers/LinearSolverParameters.cpp @@ -33,12 +33,14 @@ LinearSolverParametersInput::LinearSolverParametersInput( std::string const & na registerWrapper( viewKeyStruct::solverTypeString, &m_parameters.solverType )-> setApplyDefaultValue( m_parameters.solverType )-> setInputFlag( InputFlags::OPTIONAL )-> - setDescription( "Linear solver type. Available options are:\n* " + EnumStrings< LinearSolverParameters::SolverType >::concat( "\n* " ) ); + setDescription( "Linear solver type. Available options are:\n* " + + EnumStrings< LinearSolverParameters::SolverType >::concat( "\n* " ) ); registerWrapper( viewKeyStruct::preconditionerTypeString, &m_parameters.preconditionerType )-> setApplyDefaultValue( m_parameters.preconditionerType )-> setInputFlag( InputFlags::OPTIONAL )-> - setDescription( "Preconditioner type. Available options are:\n* " + EnumStrings< LinearSolverParameters::PreconditionerType >::concat( "\n* " ) ); + setDescription( "Preconditioner type. Available options are:\n* " + + EnumStrings< LinearSolverParameters::PreconditionerType >::concat( "\n* " ) ); registerWrapper( viewKeyStruct::stopIfErrorString, &m_parameters.stopIfError )-> setApplyDefaultValue( m_parameters.stopIfError )-> @@ -117,14 +119,14 @@ LinearSolverParametersInput::LinearSolverParametersInput( std::string const & na registerWrapper( viewKeyStruct::amgSmootherString, &m_parameters.amg.smootherType )-> setApplyDefaultValue( m_parameters.amg.smootherType )-> setInputFlag( InputFlags::OPTIONAL )-> - setDescription( "AMG smoother type\n" - "Available options are: jacobi, blockJacobi, gaussSeidel, blockGaussSeidel, chebyshev, icc, ilu, ilut" ); + setDescription( "AMG smoother type. Valid options (not all may be supported by linear algebra package):\n* " + + EnumStrings< LinearSolverParameters::PreconditionerType >::concat( "\n* " ) ); registerWrapper( viewKeyStruct::amgCoarseString, &m_parameters.amg.coarseType )-> setApplyDefaultValue( m_parameters.amg.coarseType )-> setInputFlag( InputFlags::OPTIONAL )-> - setDescription( "AMG coarsest level solver/smoother type\n" - "Available options are: jacobi, gaussSeidel, blockGaussSeidel, chebyshev, direct" ); + setDescription( "AMG coarsest level solver/smoother type. Valid options (not all may be supported by linear algebra package):\n* " + + EnumStrings< LinearSolverParameters::PreconditionerType >::concat( "\n* " ) ); registerWrapper( viewKeyStruct::amgThresholdString, &m_parameters.amg.threshold )-> setApplyDefaultValue( m_parameters.amg.threshold )-> diff --git a/src/coreComponents/physicsSolvers/SolverBase.cpp b/src/coreComponents/physicsSolvers/SolverBase.cpp index 97b2acb1542..516c1f2690a 100644 --- a/src/coreComponents/physicsSolvers/SolverBase.cpp +++ b/src/coreComponents/physicsSolvers/SolverBase.cpp @@ -241,8 +241,8 @@ real64 SolverBase::LinearImplicitStep( real64 const & time_n, // Compose parallel LA matrix/rhs out of local LA matrix/rhs m_matrix.create( m_localMatrix.toViewConst(), MPI_COMM_GEOSX ); - m_rhs.create( m_localRhs.toViewConst(), MPI_COMM_GEOSX ); - m_solution.createWithLocalSize( m_matrix.numLocalCols(), MPI_COMM_GEOSX ); + m_rhs.createWithLocalValues( m_localRhs, MPI_COMM_GEOSX ); + m_solution.createWithLocalValues( m_localSolution, MPI_COMM_GEOSX ); // Output the linear system matrix/rhs for debugging purposes DebugOutputSystem( 0.0, 0, 0, m_matrix, m_rhs ); @@ -536,8 +536,8 @@ real64 SolverBase::NonlinearImplicitStep( real64 const & time_n, // Compose parallel LA matrix/rhs out of local LA matrix/rhs m_matrix.create( m_localMatrix.toViewConst(), MPI_COMM_GEOSX ); - m_rhs.create( m_localRhs.toViewConst(), MPI_COMM_GEOSX ); - m_solution.createWithLocalSize( m_matrix.numLocalCols(), MPI_COMM_GEOSX ); + m_rhs.createWithLocalValues( m_localRhs, MPI_COMM_GEOSX ); + m_solution.createWithLocalValues( m_localSolution, MPI_COMM_GEOSX ); // Output the linear system matrix/rhs for debugging purposes DebugOutputSystem( time_n, cycleNumber, newtonIter, m_matrix, m_rhs ); @@ -641,12 +641,14 @@ void SolverBase::SetupSystem( DomainPartition & domain, dofManager.setSparsityPattern( pattern ); localMatrix.assimilate< parallelDevicePolicy<> >( std::move( pattern ) ); } + localMatrix.setName( this->getName() + "/localMatrix" ); + localRhs.move( LvArray::MemorySpace::CPU, false ); localRhs.resize( numLocalRows ); - localSolution.resize( numLocalRows ); - - localMatrix.setName( this->getName() + "/localMatrix" ); localRhs.setName( this->getName() + "/localRhs" ); + + localSolution.move( LvArray::MemorySpace::CPU, false ); + localSolution.resize( numLocalRows ); localSolution.setName( this->getName() + "/localSolution" ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/ProppantTransport.cpp b/src/coreComponents/physicsSolvers/fluidFlow/ProppantTransport.cpp index fe8eb15250b..f397e0361d0 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/ProppantTransport.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/ProppantTransport.cpp @@ -452,6 +452,7 @@ void ProppantTransport::PostStepUpdate( real64 const & time_n, UpdateProppantMobility( subRegion ); } ); + real64 const maxProppantConcentration = m_maxProppantConcentration; forTargetSubRegions( mesh, [&]( localIndex const, ElementSubRegionBase & subRegion ) { @@ -461,10 +462,10 @@ void ProppantTransport::PostStepUpdate( real64 const & time_n, forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOSX_HOST_DEVICE ( localIndex const ei ) { - if( proppantConc[ei] >= m_maxProppantConcentration || packVf[ei] >= 1.0 ) + if( proppantConc[ei] >= maxProppantConcentration || packVf[ei] >= 1.0 ) { packVf[ei] = 1.0; - proppantConc[ei] = m_maxProppantConcentration; + proppantConc[ei] = maxProppantConcentration; } } ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp index 4e14aa8bc97..fec53c45502 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp @@ -846,9 +846,9 @@ void HydrofractureSolver::ApplyBoundaryConditions( real64 const time, m_solidSolver->getSystemMatrix().create( m_solidSolver->getLocalMatrix().toViewConst(), MPI_COMM_GEOSX ); - m_solidSolver->getSystemRhs().create( m_solidSolver->getLocalRhs().toViewConst(), MPI_COMM_GEOSX ); + m_solidSolver->getSystemRhs().createWithLocalValues( m_solidSolver->getLocalRhs(), MPI_COMM_GEOSX ); m_flowSolver->getSystemMatrix().create( m_flowSolver->getLocalMatrix().toViewConst(), MPI_COMM_GEOSX ); - m_flowSolver->getSystemRhs().create( m_flowSolver->getLocalRhs().toViewConst(), MPI_COMM_GEOSX ); + m_flowSolver->getSystemRhs().createWithLocalValues( m_flowSolver->getLocalRhs(), MPI_COMM_GEOSX ); // GEOSX_LOG_RANK_0("***********************************************************"); // GEOSX_LOG_RANK_0("matrix00"); @@ -1290,14 +1290,14 @@ void HydrofractureSolver::SolveSystem( DofManager const & GEOSX_UNUSED_PARAM( do Epetra_FEVector * p_rhs[2]; Epetra_FEVector * p_solution[2]; - m_solidSolver->getSystemRhs().create( m_solidSolver->getLocalRhs(), MPI_COMM_GEOSX ); - m_flowSolver->getSystemRhs().create( m_flowSolver->getLocalRhs(), MPI_COMM_GEOSX ); + m_solidSolver->getSystemRhs().createWithLocalValues( m_solidSolver->getLocalRhs(), MPI_COMM_GEOSX ); + m_flowSolver->getSystemRhs().createWithLocalValues( m_flowSolver->getLocalRhs(), MPI_COMM_GEOSX ); p_rhs[0] = &m_solidSolver->getSystemRhs().unwrapped(); p_rhs[1] = &m_flowSolver->getSystemRhs().unwrapped(); - m_solidSolver->getSystemSolution().create( m_solidSolver->getLocalSolution(), MPI_COMM_GEOSX ); - m_flowSolver->getSystemSolution().create( m_flowSolver->getLocalSolution(), MPI_COMM_GEOSX ); + m_solidSolver->getSystemSolution().createWithLocalValues( m_solidSolver->getLocalSolution(), MPI_COMM_GEOSX ); + m_flowSolver->getSystemSolution().createWithLocalValues( m_flowSolver->getLocalSolution(), MPI_COMM_GEOSX ); p_solution[0] = &m_solidSolver->getSystemSolution().unwrapped(); p_solution[1] = &m_flowSolver->getSystemSolution().unwrapped(); diff --git a/src/coreComponents/physicsSolvers/multiphysics/LagrangianContactSolver.cpp b/src/coreComponents/physicsSolvers/multiphysics/LagrangianContactSolver.cpp index 670892c4094..57c54687c60 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/LagrangianContactSolver.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/LagrangianContactSolver.cpp @@ -630,8 +630,8 @@ real64 LagrangianContactSolver::NonlinearImplicitStep( real64 const & time_n, // Compose parallel LA matrix/rhs out of local LA matrix/rhs m_matrix.create( m_localMatrix.toViewConst(), MPI_COMM_GEOSX ); - m_rhs.create( m_localRhs.toViewConst(), MPI_COMM_GEOSX ); - m_solution.createWithLocalSize( m_matrix.numLocalCols(), MPI_COMM_GEOSX ); + m_rhs.createWithLocalValues( m_localRhs, MPI_COMM_GEOSX ); + m_solution.createWithLocalValues( m_localSolution, MPI_COMM_GEOSX ); // Output the linear system matrix/rhs for debugging purposes DebugOutputSystem( time_n, cycleNumber, newtonIter, m_matrix, m_rhs ); diff --git a/src/docs/doxygen/GeosxConfig.hpp b/src/docs/doxygen/GeosxConfig.hpp index b4840cb2401..ea9975c5a38 100644 --- a/src/docs/doxygen/GeosxConfig.hpp +++ b/src/docs/doxygen/GeosxConfig.hpp @@ -75,6 +75,8 @@ #define GEOSX_LA_INTERFACE Trilinos /// Macro defined when Trilinos interface is selected #define GEOSX_LA_INTERFACE_TRILINOS +/// Macro defined when Trilinos interface is selected +/* #undef GEOSX_LA_INTERFACE_TRILINOSTPETRA */ /// Macro defined when Hypre interface is selected /* #undef GEOSX_LA_INTERFACE_HYPRE */ /// Macro defined when PETSc interface is selected