From 9004d434d6d9cd64ddc309991dab0e43df396d06 Mon Sep 17 00:00:00 2001 From: Sergey Klevtsov Date: Fri, 25 Feb 2022 21:20:56 -0800 Subject: [PATCH] Linear algebra changes --- src/cmake/GeosxConfig.cmake | 1 + .../codingUtilities/UnitTestUtilities.hpp | 79 +++--- src/coreComponents/common/DataTypes.hpp | 8 +- src/coreComponents/common/GeosxConfig.hpp.in | 6 + src/coreComponents/common/MpiWrapper.hpp | 69 ++++++ src/coreComponents/common/TimingMacros.hpp | 6 +- .../linearAlgebra/CMakeLists.txt | 5 +- .../linearAlgebra/DofManager.cpp | 113 +++++---- .../linearAlgebra/DofManager.hpp | 47 ++-- .../linearAlgebra/common/LinearOperator.hpp | 38 ++- .../common/PreconditionerBase.hpp | 2 - .../linearAlgebra/interfaces/MatrixBase.hpp | 118 +++++++-- .../interfaces/direct/SuiteSparse.cpp | 69 +++--- .../interfaces/direct/SuperLUDist.cpp | 26 +- .../interfaces/hypre/HypreExport.cpp | 9 +- .../interfaces/hypre/HypreKernels.hpp | 190 ++++++++++----- .../interfaces/hypre/HypreMatrix.cpp | 156 +++++++++--- .../interfaces/hypre/HypreMatrix.hpp | 20 +- .../interfaces/hypre/HyprePreconditioner.cpp | 36 ++- .../interfaces/hypre/HypreUtils.cpp | 172 ++++++++++++- .../interfaces/hypre/HypreUtils.hpp | 62 ++++- .../interfaces/hypre/HypreVector.cpp | 2 +- .../interfaces/petsc/PetscMatrix.cpp | 173 ++++++++++--- .../interfaces/petsc/PetscMatrix.hpp | 13 +- .../interfaces/petsc/PetscPreconditioner.cpp | 7 +- .../interfaces/trilinos/EpetraMatrix.cpp | 230 +++++++++++++++--- .../interfaces/trilinos/EpetraMatrix.hpp | 16 ++ .../trilinos/TrilinosPreconditioner.cpp | 43 +++- .../linearAlgebra/solvers/BicgstabSolver.cpp | 10 +- .../linearAlgebra/solvers/BicgstabSolver.hpp | 6 +- .../solvers/BlockPreconditioner.cpp | 129 ++++++---- .../solvers/BlockPreconditioner.hpp | 109 ++++----- .../linearAlgebra/solvers/CgSolver.cpp | 4 +- .../linearAlgebra/solvers/CgSolver.hpp | 8 +- .../linearAlgebra/solvers/GmresSolver.cpp | 33 ++- .../linearAlgebra/solvers/GmresSolver.hpp | 11 +- .../linearAlgebra/solvers/KrylovSolver.cpp | 21 +- .../linearAlgebra/solvers/KrylovSolver.hpp | 18 +- .../linearAlgebra/solvers/KrylovUtils.hpp | 2 +- .../solvers/PreconditionerIdentity.hpp | 13 +- .../solvers/PreconditionerNull.hpp | 62 +++++ .../solvers/RichardsonSolver.cpp | 100 ++++++++ .../solvers/RichardsonSolver.hpp | 102 ++++++++ .../unitTests/testExternalSolvers.cpp | 2 +- .../unitTests/testKrylovSolvers.cpp | 27 +- .../testLinearSolverParametersEnums.cpp | 9 +- .../linearAlgebra/utilities/BlockOperator.hpp | 19 -- .../utilities/BlockOperatorView.hpp | 5 - .../utilities/InverseNormalOperator.hpp | 2 +- .../utilities/LAIHelperFunctions.hpp | 159 ++++-------- .../utilities/LinearSolverParameters.hpp | 115 ++++++++- .../utilities/NormalOperator.hpp | 19 +- .../utilities/TransposeOperator.hpp | 33 +-- src/coreComponents/mesh/CMakeLists.txt | 10 +- src/coreComponents/mesh/DomainPartition.hpp | 15 ++ .../mesh/ElementRegionManager.cpp | 3 +- .../physicsSolvers/LinearSolverParameters.cpp | 121 ++++++++- .../physicsSolvers/LinearSolverParameters.hpp | 28 ++- .../physicsSolvers/SolverBase.cpp | 47 ++-- .../physicsSolvers/SolverBase.hpp | 8 + .../contact/LagrangianContactSolver.cpp | 46 ++-- .../contact/LagrangianContactSolver.hpp | 5 +- .../fluidFlow/SinglePhaseBase.cpp | 6 - .../fluidFlow/SinglePhaseFVM.cpp | 23 +- .../fluidFlow/SinglePhaseFVM.hpp | 4 + .../multiphysics/SinglePhasePoromechanics.cpp | 121 +++++++-- .../multiphysics/SinglePhasePoromechanics.hpp | 22 +- .../SolidMechanicsLagrangianFEM.cpp | 59 +++-- .../SolidMechanicsLagrangianFEM.hpp | 30 +-- src/coreComponents/schema/docs/Block.rst | 11 + .../schema/docs/Block_other.rst | 9 + .../schema/docs/Hydrofracture.rst | 1 + .../schema/docs/LinearSolverParameters.rst | 20 +- .../docs/LinearSolverParameters_other.rst | 10 +- .../schema/docs/SinglePhasePoromechanics.rst | 1 + ...glePhasePoromechanicsEmbeddedFractures.rst | 1 + src/coreComponents/schema/schema.xsd | 64 ++++- src/coreComponents/schema/schema.xsd.other | 7 +- src/docs/doxygen/GeosxConfig.hpp | 6 + src/docs/sphinx/CompleteXMLSchema.rst | 14 ++ 80 files changed, 2486 insertions(+), 940 deletions(-) create mode 100644 src/coreComponents/linearAlgebra/solvers/PreconditionerNull.hpp create mode 100644 src/coreComponents/linearAlgebra/solvers/RichardsonSolver.cpp create mode 100644 src/coreComponents/linearAlgebra/solvers/RichardsonSolver.hpp create mode 100644 src/coreComponents/schema/docs/Block.rst create mode 100644 src/coreComponents/schema/docs/Block_other.rst diff --git a/src/cmake/GeosxConfig.cmake b/src/cmake/GeosxConfig.cmake index 124b5da4b8d..365c44eaffe 100644 --- a/src/cmake/GeosxConfig.cmake +++ b/src/cmake/GeosxConfig.cmake @@ -78,6 +78,7 @@ function( make_full_config_file set( GEOSX_LA_INTERFACE_HYPRE ON ) set( GEOSX_LA_INTERFACE_TRILINOS OFF ) set( GEOSX_LA_INTERFACE_PETSC OFF ) + set( Python3_VERSION "3.10.8" ) configure_file( ${CMAKE_SOURCE_DIR}/coreComponents/common/GeosxConfig.hpp.in ${CMAKE_SOURCE_DIR}/docs/doxygen/GeosxConfig.hpp ) diff --git a/src/coreComponents/codingUtilities/UnitTestUtilities.hpp b/src/coreComponents/codingUtilities/UnitTestUtilities.hpp index cc5280b15a7..95967a437cc 100644 --- a/src/coreComponents/codingUtilities/UnitTestUtilities.hpp +++ b/src/coreComponents/codingUtilities/UnitTestUtilities.hpp @@ -121,13 +121,11 @@ void checkRelativeError( real64 const v1, real64 const v2, real64 const relTol, EXPECT_PRED_FORMAT4( checkRelativeErrorFormat, v1, v2, relTol, absTol ); } -template< typename ROW_INDEX, typename COL_INDEX, typename VALUE > -void compareMatrixRow( ROW_INDEX const rowNumber, VALUE const relTol, VALUE const absTol, +template< typename COL_INDEX, typename VALUE > +void compareMatrixRow( VALUE const relTol, VALUE const absTol, localIndex const length1, COL_INDEX const * const indices1, VALUE const * const values1, localIndex const length2, COL_INDEX const * const indices2, VALUE const * const values2 ) { - SCOPED_TRACE( "Row " + std::to_string( rowNumber )); - EXPECT_EQ( length1, length2 ); for( localIndex j1 = 0, j2 = 0; j1 < length1 && j2 < length2; ++j1, ++j2 ) @@ -152,17 +150,31 @@ void compareMatrixRow( ROW_INDEX const rowNumber, VALUE const relTol, VALUE cons } } -template< typename ROW_INDEX, typename COL_INDEX, typename VALUE > -void compareMatrixRow( ROW_INDEX const rowNumber, VALUE const relTol, VALUE const absTol, - arraySlice1d< COL_INDEX const > indices1, arraySlice1d< VALUE const > values1, - arraySlice1d< COL_INDEX const > indices2, arraySlice1d< VALUE const > values2 ) +template< typename T, typename COL_INDEX > +void compareLocalMatrices( CRSMatrixView< T const, COL_INDEX const > const & matrix1, + CRSMatrixView< T const, COL_INDEX const > const & matrix2, + real64 const relTol = DEFAULT_REL_TOL, + real64 const absTol = DEFAULT_ABS_TOL, + globalIndex const rowOffset = 0 ) { - ASSERT_EQ( indices1.size(), values1.size() ); - ASSERT_EQ( indices2.size(), values2.size() ); + ASSERT_EQ( matrix1.numRows(), matrix2.numRows() ); + ASSERT_EQ( matrix1.numColumns(), matrix2.numColumns() ); + + matrix1.move( hostMemorySpace, false ); + matrix2.move( hostMemorySpace, false ); - compareMatrixRow( rowNumber, relTol, absTol, - indices1.size(), indices1.dataIfContiguous(), values1.dataIfContiguous(), - indices2.size(), indices2.dataIfContiguous(), values2.dataIfContiguous() ); + // check the accuracy across local rows + for( localIndex i = 0; i < matrix1.numRows(); ++i ) + { + SCOPED_TRACE( GEOS_FMT( "Row {}", i + rowOffset ) ); + compareMatrixRow( relTol, absTol, + matrix1.numNonZeros( i ), + matrix1.getColumns( i ).dataIfContiguous(), + matrix1.getEntries( i ).dataIfContiguous(), + matrix2.numNonZeros( i ), + matrix2.getColumns( i ).dataIfContiguous(), + matrix2.getEntries( i ).dataIfContiguous() ); + } } template< typename MATRIX > @@ -177,45 +189,10 @@ void compareMatrices( MATRIX const & matrix1, ASSERT_EQ( matrix1.numLocalRows(), matrix2.numLocalRows() ); ASSERT_EQ( matrix1.numLocalCols(), matrix2.numLocalCols() ); - array1d< globalIndex > indices1, indices2; - array1d< real64 > values1, values2; - - // check the accuracy across local rows - for( globalIndex i = matrix1.ilower(); i < matrix1.iupper(); ++i ) - { - indices1.resize( matrix1.rowLength( i ) ); - values1.resize( matrix1.rowLength( i ) ); - matrix1.getRowCopy( i, indices1, values1 ); - - indices2.resize( matrix2.rowLength( i ) ); - values2.resize( matrix2.rowLength( i ) ); - matrix2.getRowCopy( i, indices2, values2 ); + CRSMatrix< real64, globalIndex > const mat1 = matrix1.extract(); + CRSMatrix< real64, globalIndex > const mat2 = matrix2.extract(); - compareMatrixRow( i, relTol, absTol, - indices1.size(), indices1.data(), values1.data(), - indices2.size(), indices2.data(), values2.data() ); - } -} - -template< typename T, typename COL_INDEX > -void compareLocalMatrices( CRSMatrixView< T const, COL_INDEX const > const & matrix1, - CRSMatrixView< T const, COL_INDEX const > const & matrix2, - real64 const relTol = DEFAULT_REL_TOL, - real64 const absTol = DEFAULT_ABS_TOL ) -{ - ASSERT_EQ( matrix1.numRows(), matrix2.numRows() ); - ASSERT_EQ( matrix1.numColumns(), matrix2.numColumns() ); - - matrix1.move( hostMemorySpace, false ); - matrix2.move( hostMemorySpace, false ); - - // check the accuracy across local rows - for( localIndex i = 0; i < matrix1.numRows(); ++i ) - { - compareMatrixRow( i, relTol, absTol, - matrix1.getColumns( i ), matrix1.getEntries( i ), - matrix2.getColumns( i ), matrix2.getEntries( i ) ); - } + compareLocalMatrices( mat1.toViewConst(), mat2.toViewConst(), relTol, absTol, matrix1.ilower() ); } } // namespace testing diff --git a/src/coreComponents/common/DataTypes.hpp b/src/coreComponents/common/DataTypes.hpp index 6c81b4d58e7..9ac68ffba44 100644 --- a/src/coreComponents/common/DataTypes.hpp +++ b/src/coreComponents/common/DataTypes.hpp @@ -342,12 +342,12 @@ template< typename COL_INDEX, typename INDEX_TYPE=localIndex > using SparsityPatternView = LvArray::SparsityPatternView< COL_INDEX, INDEX_TYPE const, LvArray::ChaiBuffer >; /// Alias for CRS Matrix class. -template< typename T, typename COL_INDEX=globalIndex > -using CRSMatrix = LvArray::CRSMatrix< T, COL_INDEX, localIndex, LvArray::ChaiBuffer >; +template< typename T, typename COL_INDEX=globalIndex, typename INDEX_TYPE=localIndex > +using CRSMatrix = LvArray::CRSMatrix< T, COL_INDEX, INDEX_TYPE, LvArray::ChaiBuffer >; /// Alias for CRS Matrix View. -template< typename T, typename COL_INDEX=globalIndex > -using CRSMatrixView = LvArray::CRSMatrixView< T, COL_INDEX, localIndex const, LvArray::ChaiBuffer >; +template< typename T, typename COL_INDEX=globalIndex, typename INDEX_TYPE=localIndex > +using CRSMatrixView = LvArray::CRSMatrixView< T, COL_INDEX, INDEX_TYPE const, LvArray::ChaiBuffer >; ///@} diff --git a/src/coreComponents/common/GeosxConfig.hpp.in b/src/coreComponents/common/GeosxConfig.hpp.in index 2e0bad1a7f9..f582f045e93 100644 --- a/src/coreComponents/common/GeosxConfig.hpp.in +++ b/src/coreComponents/common/GeosxConfig.hpp.in @@ -92,6 +92,12 @@ /// Enables use of PETSc library (CMake option ENABLE_PETSC) #cmakedefine GEOSX_USE_PETSC +/// Enables use of METIS library (CMake option ENABLE_METIS) +#cmakedefine GEOSX_USE_METIS + +/// Enables use of ParMETIS library (CMake option ENABLE_PARMETIS) +#cmakedefine GEOSX_USE_PARMETIS + /// Enables use of Scotch library (CMake option ENABLE_SCOTCH) #cmakedefine GEOSX_USE_SCOTCH diff --git a/src/coreComponents/common/MpiWrapper.hpp b/src/coreComponents/common/MpiWrapper.hpp index 5dde9e71200..f11fb8b5df6 100644 --- a/src/coreComponents/common/MpiWrapper.hpp +++ b/src/coreComponents/common/MpiWrapper.hpp @@ -22,6 +22,8 @@ #include "common/Span.hpp" #include "mesh/ElementType.hpp" +#include + #if defined(GEOSX_USE_MPI) #include #define MPI_PARAM( x ) x @@ -343,6 +345,11 @@ struct MpiWrapper array1d< T > & recvbuf, MPI_Comm comm = MPI_COMM_GEOSX ); + template< typename T > + static int allGatherv( arrayView1d< T const > const & sendbuf, + array1d< T > & recvbuf, + MPI_Comm comm = MPI_COMM_GEOSX ); + /** * @brief Strongly typed wrapper around MPI_Allreduce. * @param[in] sendbuf The pointer to the sending buffer. @@ -465,6 +472,22 @@ struct MpiWrapper MPI_Comm comm, MPI_Request * request ); + /** + * @brief Strongly typed wrapper around MPI_Send() + * @param[in] buf The pointer to the buffer that contains the data to be sent. + * @param[in] count The number of elements in \p buf. + * @param[in] dest The rank of the destination process within \p comm. + * @param[in] tag The message tag that is be used to distinguish different types of messages. + * @param[in] comm The handle to the MPI_Comm. + * @return + */ + template< typename T > + static int send( T const * const buf, + int count, + int dest, + int tag, + MPI_Comm comm ); + /** * @brief Strongly typed wrapper around MPI_Isend() * @param[in] buf The pointer to the buffer that contains the data to be sent. @@ -719,6 +742,38 @@ int MpiWrapper::allGather( arrayView1d< T const > const & sendValues, #endif } +template< typename T > +int MpiWrapper::allGatherv( arrayView1d< T const > const & sendValues, + array1d< T > & allValues, + MPI_Comm MPI_PARAM( comm ) ) +{ + int const sendSize = LvArray::integerConversion< int >( sendValues.size() ); +#ifdef GEOSX_USE_MPI + int const mpiSize = commSize( comm ); + array1d< int > counts; + allGather( sendSize, counts, comm ); + array1d< int > displs( mpiSize + 1 ); + std::partial_sum( counts.begin(), counts.end(), displs.begin() + 1 ); + allValues.resize( displs.back() ); + return MPI_Allgatherv( sendValues.data(), + sendSize, + internal::getMpiType< T >(), + allValues.data(), + counts.data(), + displs.data(), + internal::getMpiType< T >(), + comm ); + +#else + allValues.resize( sendSize ); + for( localIndex a=0; a int MpiWrapper::allReduce( T const * const sendbuf, T * const recvbuf, @@ -934,6 +989,20 @@ int MpiWrapper::iSend( arrayView1d< T > const & buf, #endif } +template< typename T > +int MpiWrapper::send( T const * const buf, + int count, + int dest, + int tag, + MPI_Comm comm ) +{ +#ifdef GEOSX_USE_MPI + return MPI_Send( buf, count, internal::getMpiType< T >(), dest, tag, comm ); +#else + GEOS_ERROR( "Not implemented without MPI" ); +#endif +} + template< typename T > int MpiWrapper::iSend( T const * const buf, int count, diff --git a/src/coreComponents/common/TimingMacros.hpp b/src/coreComponents/common/TimingMacros.hpp index 297aa49f8f1..d1f162bd6b4 100644 --- a/src/coreComponents/common/TimingMacros.hpp +++ b/src/coreComponents/common/TimingMacros.hpp @@ -34,7 +34,7 @@ namespace timingHelpers { std::string input(prettyFunction); std::string::size_type const end = input.find_first_of( '(' ); - std::string::size_type const beg = input.find_last_of( ' ', end)+1; + std::string::size_type const beg = input.find_last_of( ' ', end ) + 1; return input.substr( beg, end-beg ); } } @@ -95,10 +95,10 @@ namespace timingHelpers #include /// Mark a function or scope for timing with a given name -#define GEOS_CALIPER_MARK_SCOPE(name) cali::Function __cali_ann##__LINE__(STRINGIZE_NX(name)) +#define GEOS_CALIPER_MARK_SCOPE(name) cali::Function GEOS_CONCAT(_cali_ann, __LINE__)(STRINGIZE_NX(name)) /// Mark a function for timing using a compiler-provided name -#define GEOS_CALIPER_MARK_FUNCTION cali::Function __cali_ann##__func__(timingHelpers::stripPF(__PRETTY_FUNCTION__).c_str()) +#define GEOS_CALIPER_MARK_FUNCTION cali::Function _cali_ann_func(timingHelpers::stripPF(__PRETTY_FUNCTION__).c_str()) /// Mark the beginning of timed statement group #define GEOS_CALIPER_MARK_BEGIN(name) CALI_MARK_BEGIN(STRINGIZE(name)) diff --git a/src/coreComponents/linearAlgebra/CMakeLists.txt b/src/coreComponents/linearAlgebra/CMakeLists.txt index 271cb233b3c..1ff242e1a8f 100644 --- a/src/coreComponents/linearAlgebra/CMakeLists.txt +++ b/src/coreComponents/linearAlgebra/CMakeLists.txt @@ -21,6 +21,8 @@ set( linearAlgebra_headers solvers/PreconditionerBlockJacobi.hpp solvers/PreconditionerIdentity.hpp solvers/PreconditionerJacobi.hpp + solvers/PreconditionerNull.hpp + solvers/RichardsonSolver.hpp solvers/SeparateComponentPreconditioner.hpp utilities/Arnoldi.hpp utilities/BlockOperator.hpp @@ -49,8 +51,9 @@ set( linearAlgebra_sources solvers/CgSolver.cpp solvers/GmresSolver.cpp solvers/KrylovSolver.cpp + solvers/RichardsonSolver.cpp solvers/SeparateComponentPreconditioner.cpp - utilities/ReverseCutHillMcKeeOrdering.cpp + utilities/ReverseCutHillMcKeeOrdering.cpp ) set( dependencyList ${parallelDeps} mesh denseLinearAlgebra ) diff --git a/src/coreComponents/linearAlgebra/DofManager.cpp b/src/coreComponents/linearAlgebra/DofManager.cpp index 12380d90a59..1058c76f469 100644 --- a/src/coreComponents/linearAlgebra/DofManager.cpp +++ b/src/coreComponents/linearAlgebra/DofManager.cpp @@ -63,11 +63,11 @@ void DofManager::setDomain( DomainPartition & domain ) m_domain = &domain; } -localIndex DofManager::getFieldIndex( string const & name ) const +integer DofManager::getFieldIndex( string const & name ) const { auto const it = std::find_if( m_fields.begin(), m_fields.end(), [&]( FieldDescription const & f ) { return f.name == name; } ); - GEOS_ASSERT_MSG( it != m_fields.end(), "DofManager: field does not exist: " << name ); + GEOS_ERROR_IF( it == m_fields.end(), "DofManager: field does not exist: " << name ); return std::distance( m_fields.begin(), it ); } @@ -145,6 +145,11 @@ globalIndex DofManager::globalOffset( string const & fieldName ) const return m_fields[getFieldIndex( fieldName )].globalOffset; } +Span< DofManager::FieldSupport const > DofManager::support( string const & fieldName ) const +{ + return m_fields[getFieldIndex( fieldName )].support; +} + namespace { @@ -330,7 +335,7 @@ void DofManager::removeIndexArray( FieldDescription const & field ) } ); } -void DofManager::computeFieldDimensions( localIndex const fieldIndex ) +void DofManager::computeFieldDimensions( integer const fieldIndex ) { FieldDescription & field = m_fields[fieldIndex]; @@ -432,7 +437,7 @@ void DofManager::addField( string const & fieldName, { GEOS_ASSERT_MSG( m_domain != nullptr, "Domain has not been set" ); GEOS_ERROR_IF( m_reordered, "Cannot add fields after reorderByRank() has been called." ); - GEOS_ERROR_IF_GT_MSG( components, MAX_COMP, "Number of components limit exceeded" ); + GEOS_ERROR_IF_GT_MSG( components, maxNumComp, "Number of components limit exceeded" ); std::vector< FieldSupport > processedSupports = processFieldSupportList( *m_domain, regions ); @@ -602,8 +607,8 @@ void DofManager::addCoupling( string const & rowFieldName, bool const symmetric ) { GEOS_ASSERT_MSG( m_domain != nullptr, "Domain has not been set" ); - localIndex const rowFieldIndex = getFieldIndex( rowFieldName ); - localIndex const colFieldIndex = getFieldIndex( colFieldName ); + integer const rowFieldIndex = getFieldIndex( rowFieldName ); + integer const colFieldIndex = getFieldIndex( colFieldName ); // Check if already defined std::vector< FieldSupport > processSupportList = processCouplingRegionList( supports, @@ -646,7 +651,7 @@ void DofManager::addCoupling( string const & rowFieldName, void DofManager::addCoupling( string const & fieldName, FluxApproximationBase const & stencils ) { - localIndex const fieldIndex = getFieldIndex( fieldName ); + integer const fieldIndex = getFieldIndex( fieldName ); FieldDescription const & field = m_fields[fieldIndex]; GEOS_ERROR_IF( field.location != FieldLocation::Elem, "Field must be supported on elements in order to use stencil sparsity" ); @@ -826,11 +831,11 @@ void makeConnLocPattern( MeshLevel const & mesh, } // namespace void DofManager::setSparsityPatternFromStencil( SparsityPatternView< globalIndex > const & pattern, - localIndex const fieldIndex ) const + integer const fieldIndex ) const { FieldDescription const & field = m_fields[fieldIndex]; CouplingDescription const & coupling = m_coupling.at( {fieldIndex, fieldIndex} ); - localIndex const numComp = field.numComponents; + integer const numComp = field.numComponents; globalIndex const rankDofOffset = rankOffset(); CompMask const & globallyCoupledComps = field.globallyCoupledComponents; @@ -872,7 +877,7 @@ void DofManager::setSparsityPatternFromStencil( SparsityPatternView< globalIndex colDofIndices.resize( stencilSize * numComp ); for( localIndex i = 0; i < stencilSize; ++i ) { - for( localIndex c = 0; c < numComp; ++c ) + for( integer c = 0; c < numComp; ++c ) { colDofIndices[i * numComp + c] = dofNumber[seri( iconn, i )][sesri( iconn, i )][sei( iconn, i )] + c; } @@ -901,11 +906,11 @@ void DofManager::setSparsityPatternFromStencil( SparsityPatternView< globalIndex forMeshLocation< FieldLocation::Elem, false, parallelHostPolicy >( mesh, regions, [=]( auto const & elemIdx ) { globalIndex const elemDof = dofNumberView[elemIdx[0]][elemIdx[1]][elemIdx[2]]; - for( localIndex c = 0; c < numComp; ++c ) + for( integer c = 0; c < numComp; ++c ) { colDofIndices[c] = elemDof + c; } - for( localIndex c = 0; c < numComp; ++c ) + for( integer c = 0; c < numComp; ++c ) { pattern.insertNonZeros( elemDof - rankDofOffset + c, colDofIndices.begin(), colDofIndices.end() ); } @@ -914,8 +919,8 @@ void DofManager::setSparsityPatternFromStencil( SparsityPatternView< globalIndex } void DofManager::setSparsityPatternOneBlock( SparsityPatternView< globalIndex > const & pattern, - localIndex const rowFieldIndex, - localIndex const colFieldIndex ) const + integer const rowFieldIndex, + integer const colFieldIndex ) const { GEOS_ASSERT( rowFieldIndex >= 0 ); GEOS_ASSERT( colFieldIndex >= 0 ); @@ -1041,14 +1046,14 @@ void DofManager::setSparsityPatternOneBlock( SparsityPatternView< globalIndex > } void DofManager::countRowLengthsFromStencil( arrayView1d< localIndex > const & rowLengths, - localIndex const fieldIndex ) const + integer const fieldIndex ) const { FieldDescription const & field = m_fields[fieldIndex]; CouplingDescription const & coupling = m_coupling.at( {fieldIndex, fieldIndex} ); GEOS_ASSERT( coupling.connector == Connector::Stencil ); GEOS_ASSERT( coupling.stencils != nullptr ); - localIndex const numComp = field.numComponents; + integer const numComp = field.numComponents; globalIndex const rankDofOffset = rankOffset(); CompMask const & globallyCoupledComps = field.globallyCoupledComponents; @@ -1096,7 +1101,7 @@ void DofManager::countRowLengthsFromStencil( arrayView1d< localIndex > const & r if( ghostRank[ei] < 0 ) { localIndex const localDofNumber = dofNumber[ei] - rankDofOffset; - for( localIndex c = 0; c < numComp; ++c ) + for( integer c = 0; c < numComp; ++c ) { rowLengths[localDofNumber + c] += numComp; } @@ -1107,8 +1112,8 @@ void DofManager::countRowLengthsFromStencil( arrayView1d< localIndex > const & r } void DofManager::countRowLengthsOneBlock( arrayView1d< localIndex > const & rowLengths, - localIndex const rowFieldIndex, - localIndex const colFieldIndex ) const + integer const rowFieldIndex, + integer const colFieldIndex ) const { GEOS_ASSERT( rowFieldIndex >= 0 ); GEOS_ASSERT( colFieldIndex >= 0 ); @@ -1236,13 +1241,13 @@ void DofManager::setSparsityPattern( SparsityPattern< globalIndex > & pattern ) GEOS_ERROR_IF( !m_reordered, "Cannot set monolithic sparsity pattern before reorderByRank() has been called." ); localIndex const numLocalRows = numLocalDofs(); - localIndex const numFields = LvArray::integerConversion< localIndex >( m_fields.size() ); + integer const numFields = LvArray::integerConversion< integer >( m_fields.size() ); // Step 1. Do a dry run of sparsity construction to get the total number of nonzeros in each row array1d< localIndex > rowSizes( numLocalRows ); - for( localIndex blockRow = 0; blockRow < numFields; ++blockRow ) + for( integer blockRow = 0; blockRow < numFields; ++blockRow ) { - for( localIndex blockCol = 0; blockCol < numFields; ++blockCol ) + for( integer blockCol = 0; blockCol < numFields; ++blockCol ) { countRowLengthsOneBlock( rowSizes, blockRow, blockCol ); } @@ -1252,9 +1257,9 @@ void DofManager::setSparsityPattern( SparsityPattern< globalIndex > & pattern ) pattern.resizeFromRowCapacities< parallelHostPolicy >( numLocalRows, numGlobalDofs(), rowSizes.data() ); // Step 3. Fill the sparsity block-by-block - for( localIndex blockRow = 0; blockRow < numFields; ++blockRow ) + for( integer blockRow = 0; blockRow < numFields; ++blockRow ) { - for( localIndex blockCol = 0; blockCol < numFields; ++blockCol ) + for( integer blockCol = 0; blockCol < numFields; ++blockCol ) { setSparsityPatternOneBlock( pattern.toView(), blockRow, blockCol ); } @@ -1578,7 +1583,7 @@ void DofManager::reorderByRank() for( FieldDescription & field : m_fields ) { // compute field dimensions (local dofs, global dofs, etc) - computeFieldDimensions( static_cast< localIndex >( getFieldIndex( field.name ) ) ); + computeFieldDimensions( static_cast< integer >( getFieldIndex( field.name ) ) ); // allocate and fill index array createIndexArray( field, permutations.at( field.name ).toViewConst() ); } @@ -1592,7 +1597,7 @@ void DofManager::reorderByRank() } // This is a map with a key that is the pair of strings specifying the - // ( MeshBody name, MeshLevel name), and a value that is another map with a + // (MeshBody name, MeshLevel name), and a value that is another map with a // key that indicates the name of the object that contains the field to be // synced, and a value that contans the name of the field to be synced. std::map< std::pair< string, string >, FieldIdentifiers > fieldsToBeSync; @@ -1661,6 +1666,7 @@ void DofManager::setupFrom( DofManager const & source, std::vector< SubComponent > const & selection ) { clear(); + m_domain = source.m_domain; for( FieldDescription const & field : source.m_fields ) { auto const it = std::find_if( selection.begin(), selection.end(), @@ -1678,8 +1684,8 @@ void DofManager::setupFrom( DofManager const & source, string const & colFieldName = source.m_fields[entry.first.second].name; if( fieldExists( rowFieldName ) && fieldExists( colFieldName ) ) { - localIndex const rowFieldIndex = getFieldIndex( rowFieldName ); - localIndex const colFieldIndex = getFieldIndex( colFieldName ); + integer const rowFieldIndex = getFieldIndex( rowFieldName ); + integer const colFieldIndex = getFieldIndex( colFieldName ); m_coupling.insert( { {rowFieldIndex, colFieldIndex}, entry.second } ); } } @@ -1695,7 +1701,6 @@ void DofManager::makeRestrictor( std::vector< SubComponent > const & selection, GEOS_ERROR_IF( !m_reordered, "Cannot make restrictors before reorderByRank() has been called." ); // 1. Populate selected fields and compute some basic dimensions - // array1d< FieldDescription > fieldsSelected( selection.size() ); std::vector< FieldDescription > fieldsSelected( selection.size() ); for( std::size_t k = 0; k < fieldsSelected.size(); ++k ) @@ -1720,7 +1725,7 @@ void DofManager::makeRestrictor( std::vector< SubComponent > const & selection, blockOffset += field.numGlobalDof; } - globalIndex globalOffset = std::accumulate( fieldsSelected.begin(), fieldsSelected.end(), globalIndex( 0 ), + globalIndex globalOffset = std::accumulate( fieldsSelected.begin(), fieldsSelected.end(), globalIndex{}, []( localIndex const n, FieldDescription const & f ) { return n + f.rankOffset; } ); @@ -1732,7 +1737,7 @@ void DofManager::makeRestrictor( std::vector< SubComponent > const & selection, // 3. Build the restrictor field by field - localIndex const numLocalDofSelected = std::accumulate( fieldsSelected.begin(), fieldsSelected.end(), localIndex( 0 ), + localIndex const numLocalDofSelected = std::accumulate( fieldsSelected.begin(), fieldsSelected.end(), localIndex{}, []( localIndex const n, FieldDescription const & f ) { return n + f.numLocalDof; } ); @@ -1744,44 +1749,44 @@ void DofManager::makeRestrictor( std::vector< SubComponent > const & selection, array1d< globalIndex > rows; array1d< globalIndex > cols; - array1d< real64 > values; + array1d< real64 > vals; for( std::size_t k = 0; k < fieldsSelected.size(); ++k ) { FieldDescription const & fieldNew = fieldsSelected[k]; FieldDescription const & fieldOld = m_fields[getFieldIndex( fieldNew.name )]; - FieldDescription const & fieldRow = transpose ? fieldOld : fieldNew; - FieldDescription const & fieldCol = transpose ? fieldNew : fieldOld; - CompMask const mask = selection[k].mask; localIndex const numLocalNodes = fieldNew.numLocalDof / fieldNew.numComponents; + rows.resize( fieldNew.numLocalDof ); + cols.resize( fieldNew.numLocalDof ); + vals.resize( fieldNew.numLocalDof ); + vals.setValues< parallelHostPolicy >( 1.0 ); - rows.resize( numLocalNodes*mask.size() ); - cols.resize( numLocalNodes*mask.size() ); - values.resize( numLocalNodes*mask.size() ); - - for( localIndex i = 0; i < numLocalNodes; ++i ) + forAll< parallelHostPolicy >( numLocalNodes, [&fieldNew, &fieldOld, mask, transpose, + rows = rows.toView(), + cols = cols.toView()]( localIndex const i ) { integer newComp = 0; for( integer const oldComp : mask ) { - globalIndex const row = fieldRow.globalOffset + i * fieldRow.numComponents + ( transpose ? oldComp : newComp ); - globalIndex const col = fieldCol.globalOffset + i * fieldCol.numComponents + ( transpose ? newComp : oldComp ); + globalIndex row = fieldNew.globalOffset + i * fieldNew.numComponents + newComp; + globalIndex col = fieldOld.globalOffset + i * fieldOld.numComponents + oldComp; + if( transpose ) + { + std::swap( row, col ); + } - rows( i*mask.size() + newComp ) = row; - cols( i*mask.size() + newComp ) = col; - values[ i*mask.size() + newComp ] = 1.0; + rows[ i * fieldNew.numComponents + newComp ] = row; + cols[ i * fieldNew.numComponents + newComp ] = col; ++newComp; } - } + } ); restrictor.insert( rows.toViewConst(), cols.toViewConst(), - values.toViewConst() ); - - + vals.toViewConst() ); } restrictor.close(); } @@ -1791,12 +1796,12 @@ void DofManager::printFieldInfo( std::ostream & os ) const { if( MpiWrapper::commRank( MPI_COMM_GEOSX ) == 0 ) { - localIndex const numFields = LvArray::integerConversion< localIndex >( m_fields.size() ); + integer const numFields = LvArray::integerConversion< integer >( m_fields.size() ); os << "Fields:" << std::endl; os << " # | " << std::setw( 20 ) << "name" << " | " << "comp" << " | " << "N global DOF" << std::endl; os << "---+----------------------+------+-------------" << std::endl; - for( localIndex i = 0; i < numFields; ++i ) + for( integer i = 0; i < numFields; ++i ) { FieldDescription const & f = m_fields[i]; os << ' ' << i << " | " @@ -1807,9 +1812,9 @@ void DofManager::printFieldInfo( std::ostream & os ) const os << "---+----------------------+------+-------------" << std::endl; os << std::endl << "Connectivity:" << std::endl; - for( localIndex i = 0; i < numFields; ++i ) + for( integer i = 0; i < numFields; ++i ) { - for( localIndex j = 0; j < numFields; ++j ) + for( integer j = 0; j < numFields; ++j ) { if( m_coupling.count( {i, j} ) == 0 ) { @@ -1860,7 +1865,7 @@ void DofManager::printFieldInfo( std::ostream & os ) const os << std::endl; if( i < numFields - 1 ) { - for( localIndex j = 0; j < numFields - 1; ++j ) + for( integer j = 0; j < numFields - 1; ++j ) { os << "---|"; } diff --git a/src/coreComponents/linearAlgebra/DofManager.hpp b/src/coreComponents/linearAlgebra/DofManager.hpp index db8b7a01c8d..b544ae22c47 100644 --- a/src/coreComponents/linearAlgebra/DofManager.hpp +++ b/src/coreComponents/linearAlgebra/DofManager.hpp @@ -20,6 +20,7 @@ #define GEOS_LINEARALGEBRA_DOFMANAGER_HPP_ #include "common/DataTypes.hpp" +#include "common/Span.hpp" #include "linearAlgebra/utilities/ComponentMask.hpp" #include "mesh/FieldIdentifiers.hpp" @@ -44,10 +45,10 @@ class DofManager public: /// Maximum number of components in a field - static int constexpr MAX_COMP = 32; + static int constexpr maxNumComp = 32; /// Type of component mask used by DofManager - using CompMask = ComponentMask< MAX_COMP >; + using CompMask = ComponentMask< maxNumComp >; /** * @brief Describes a selection of components from a DoF field. @@ -325,7 +326,7 @@ class DofManager * @brief @return number of components in a given field. * @param fieldName name of the field */ - integer numComponents( string const & fieldName = "" ) const; + integer numComponents( string const & fieldName ) const; /** * @brief @return number of dof components across all fields. @@ -345,6 +346,12 @@ class DofManager */ globalIndex globalOffset( string const & fieldName ) const; + /** + * @brief @return support of the given field (list of mesh body/levels/regions) + * @param fieldName + */ + Span< FieldSupport const > support( string const & fieldName ) const; + /** * @brief Return an array of number of components per field, sorted by field registration order. * @return array of number of components @@ -367,7 +374,7 @@ class DofManager auto it = labels.begin(); for( FieldDescription const & field : m_fields ) { - localIndex const numComp = field.numComponents; + integer const numComp = field.numComponents; localIndex const numSupp = field.numLocalDof / numComp; for( localIndex i = 0; i < numSupp; ++i, it += numComp ) { @@ -396,7 +403,7 @@ class DofManager string const & srcFieldName, string const & dstFieldName, real64 scalingFactor, - CompMask mask = CompMask( MAX_COMP, true ) ) const; + CompMask mask = CompMask( maxNumComp, true ) ) const; /** * @brief Add values from LA vectors to simulation data arrays. @@ -412,7 +419,7 @@ class DofManager string const & srcFieldName, string const & dstFieldName, SCALING_FACTOR_TYPE const & scalingFactor, - CompMask mask = CompMask( MAX_COMP, true ) ) const; + CompMask mask = CompMask( maxNumComp, true ) ) const; /** * @brief Copy values from simulation data arrays to vectors. @@ -427,7 +434,7 @@ class DofManager string const & srcFieldName, string const & dstFieldName, real64 scalingFactor, - CompMask mask = CompMask( MAX_COMP, true ) ) const; + CompMask mask = CompMask( maxNumComp, true ) ) const; /** * @brief Add values from a simulation data array to a DOF vector. @@ -442,7 +449,7 @@ class DofManager string const & srcFieldName, string const & dstFieldName, real64 scalingFactor, - CompMask mask = CompMask( MAX_COMP, true ) ) const; + CompMask mask = CompMask( maxNumComp, true ) ) const; /** * @brief Create a dof selection by filtering out excluded components @@ -526,13 +533,13 @@ class DofManager /** * @brief Get field index from string key */ - localIndex getFieldIndex( string const & name ) const; + integer getFieldIndex( string const & name ) const; /** * @brief Compute and save dof offsets the field * @param fieldIndex index of the field */ - void computeFieldDimensions( localIndex fieldIndex ); + void computeFieldDimensions( integer const fieldIndex ); /** * @brief Create index array for the field @@ -572,11 +579,11 @@ class DofManager * @param colFieldIndex index of col field (must be non-negative) */ void countRowLengthsOneBlock( arrayView1d< localIndex > const & rowLengths, - localIndex rowFieldIndex, - localIndex colFieldIndex ) const; + integer rowFieldIndex, + integer colFieldIndex ) const; void countRowLengthsFromStencil( arrayView1d< localIndex > const & rowLengths, - localIndex fieldIndex ) const; + integer fieldIndex ) const; /** * @brief Populate the sparsity pattern for a coupling block between given fields. @@ -587,15 +594,11 @@ class DofManager * This private function is used as a building block by higher-level SetSparsityPattern() */ void setSparsityPatternOneBlock( SparsityPatternView< globalIndex > const & pattern, - localIndex rowFieldIndex, - localIndex colFieldIndex ) const; + integer rowFieldIndex, + integer colFieldIndex ) const; void setSparsityPatternFromStencil( SparsityPatternView< globalIndex > const & pattern, - localIndex fieldIndex ) const; - - template< int DIMS_PER_DOF > - void setFiniteElementSparsityPattern( SparsityPattern< globalIndex > & pattern, - localIndex fieldIndex ) const; + integer fieldIndex ) const; /** * @brief Generic implementation for @ref copyVectorToField and @ref addVectorToField @@ -634,14 +637,14 @@ class DofManager /// Name of the manager (unique, for unique identification of index array keys) string m_name; - /// Pointer to corresponding MeshLevel + /// Pointer to domain that holds mesh bodies/levels DomainPartition * m_domain = nullptr; /// Array of field descriptions std::vector< FieldDescription > m_fields; /// Table of connector types within and between fields - std::map< std::pair< localIndex, localIndex >, CouplingDescription > m_coupling; + std::map< std::pair< integer, integer >, CouplingDescription > m_coupling; /// Flag indicating that DOFs have been reordered rank-wise. bool m_reordered = false; diff --git a/src/coreComponents/linearAlgebra/common/LinearOperator.hpp b/src/coreComponents/linearAlgebra/common/LinearOperator.hpp index 995ef2b6936..8dc5fc5ca99 100644 --- a/src/coreComponents/linearAlgebra/common/LinearOperator.hpp +++ b/src/coreComponents/linearAlgebra/common/LinearOperator.hpp @@ -38,12 +38,12 @@ class LinearOperator using Vector = VECTOR; /** - * @brief Constructor + * @brief Constructor. */ LinearOperator() = default; /** - * @brief Destructor + * @brief Destructor. */ virtual ~LinearOperator() = default; @@ -73,26 +73,36 @@ class LinearOperator } /** - * @brief Get the number of global rows. - * @return Number of global rows in the operator. + * @brief @return the number of global rows. */ virtual globalIndex numGlobalRows() const = 0; /** - * @brief Get the number of global columns. - * @return Number of global columns in the operator. + * @brief @return the number of global columns. */ virtual globalIndex numGlobalCols() const = 0; /** - * @brief Get the number of local rows. - * @return Number of local rows in the operator. + * @brief @return the total number of nonzero entries in the operator. + * + * @note While "non-zero entries" is not a well-defined term for the abstract concept of a linear operator, + * in practice this value is needed to track runtime/memory complexity of the operator implementation. + * The function returns the number of nonzero floating-point values stored/used when applying to a vector. + * A matrix-free operator implemented without any additional FP storage may either return 0, + * or the (approximate) number of floating-point multiply operations performed by apply(). + */ + virtual globalIndex numGlobalNonzeros() const + { + return 0; + } + + /** + * @brief @return the number of local rows. */ virtual localIndex numLocalRows() const = 0; /** - * @brief Get the number of local columns. - * @return Number of local columns in the operator. + * @brief @return the number of local columns. * * @note The use of term "local columns" refers not to physical partitioning of columns across ranks * (as e.g. matrices are partitioned by rows and typically physically store all column entries), @@ -100,6 +110,14 @@ class LinearOperator */ virtual localIndex numLocalCols() const = 0; + /** + * @brief @return the number of nonzero entries in the local portion of the operator. + */ + virtual localIndex numLocalNonzeros() const + { + return 0; + } + /** * @brief Get the MPI communicator the matrix was created with * @return MPI communicator passed in @p create...() diff --git a/src/coreComponents/linearAlgebra/common/PreconditionerBase.hpp b/src/coreComponents/linearAlgebra/common/PreconditionerBase.hpp index d8869c7eb2e..a65ba5fd1f4 100644 --- a/src/coreComponents/linearAlgebra/common/PreconditionerBase.hpp +++ b/src/coreComponents/linearAlgebra/common/PreconditionerBase.hpp @@ -34,8 +34,6 @@ class PreconditionerBase : public LinearOperator< typename LAI::ParallelVector > PreconditionerBase() = default; - virtual ~PreconditionerBase() = default; - /// Alias for base type using Base = LinearOperator< typename LAI::ParallelVector >; diff --git a/src/coreComponents/linearAlgebra/interfaces/MatrixBase.hpp b/src/coreComponents/linearAlgebra/interfaces/MatrixBase.hpp index 65a9df9b176..2c3da50decd 100644 --- a/src/coreComponents/linearAlgebra/interfaces/MatrixBase.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/MatrixBase.hpp @@ -45,7 +45,7 @@ enum class RowSumType */ enum class MatrixPatternOp { - Same, // Caller guarantees patterns of arguments are exactly the same + Equal, // Caller guarantees patterns of arguments are exactly the same Subset, // Caller guarantees pattern of second argument is a subset of the first Restrict, // Restrict pattern of second argument, ignoring any entries that don't exist in first Extend // Extend pattern of the first argument with entries from the second @@ -84,10 +84,12 @@ class MatrixBase : public virtual LinearOperator< VECTOR > /// Type alias for a compatible vector class using Vector = VECTOR; - using Base::numLocalRows; - using Base::numLocalCols; using Base::numGlobalRows; using Base::numGlobalCols; + using Base::numGlobalNonzeros; + using Base::numLocalRows; + using Base::numLocalCols; + using Base::numLocalNonzeros; /** * @name Status query methods @@ -616,6 +618,27 @@ class MatrixBase : public virtual LinearOperator< VECTOR > R.multiply( AP, dst ); } + /** + * @brief Compute the triple product dst = P1^T * this * P2 + * @param P1 the left "prolongation" matrix + * @param P2 the right "prolongation" matrix + * @param dst the resulting product matrix (will be re-created as needed) + */ + virtual void multiplyPtAP( Matrix const & P1, + Matrix const & P2, + Matrix & dst ) const + { + GEOS_LAI_ASSERT( ready() ); + GEOS_LAI_ASSERT( P1.ready() ); + GEOS_LAI_ASSERT( P2.ready() ); + GEOS_LAI_ASSERT_EQ( numLocalRows(), P1.numLocalRows() ); + GEOS_LAI_ASSERT_EQ( numLocalCols(), P2.numLocalRows() ); + + Matrix AP2; + multiply( P2, AP2 ); + P1.leftMultiplyTranspose( AP2, dst ); + } + /** * @brief Compute the triple product dst = P^T * this * P * @param P the "prolongation" matrix @@ -624,14 +647,7 @@ class MatrixBase : public virtual LinearOperator< VECTOR > virtual void multiplyPtAP( Matrix const & P, Matrix & dst ) const { - GEOS_LAI_ASSERT( ready() ); - GEOS_LAI_ASSERT( P.ready() ); - GEOS_LAI_ASSERT_EQ( numGlobalRows(), P.numGlobalRows() ); - GEOS_LAI_ASSERT_EQ( numGlobalCols(), P.numGlobalRows() ); - - Matrix AP; - multiply( P, AP ); - P.leftMultiplyTranspose( AP, dst ); + multiplyPtAP( P, P, dst ); } /** @@ -777,13 +793,24 @@ class MatrixBase : public virtual LinearOperator< VECTOR > */ ///@{ + /** + * @brief Returns the number of nonzero entries in the longest local row of the matrix. + * @return the max length of a row + * + * Not collective. + */ + virtual localIndex maxRowLengthLocal() const = 0; + /** * @brief Returns the number of nonzero entries in the longest row of the matrix. * @return the max length of a row * * Collective. */ - virtual localIndex maxRowLength() const = 0; + virtual localIndex maxRowLength() const + { + return MpiWrapper::max( maxRowLengthLocal(), this->comm() ); + } /** * @brief Get row length via global row index. @@ -799,6 +826,13 @@ class MatrixBase : public virtual LinearOperator< VECTOR > */ virtual void getRowLengths( arrayView1d< localIndex > const & lengths ) const = 0; + /** + * @brief Get the local row lengths of every local row (number of nonzeros in local columns) + * @param lengths an array view to be populated with row lengths + * @note The implementation may move the view's buffer to a different memory space. + */ + virtual void getRowLocalLengths( arrayView1d< localIndex > const & lengths ) const = 0; + /** * @brief Returns a copy of the data in row @p globalRow. * @param[in] globalRow the index of global row to extract @@ -815,6 +849,54 @@ class MatrixBase : public virtual LinearOperator< VECTOR > */ virtual void extractDiagonal( Vector & dst ) const = 0; + /** + * @brief Extract local part of the matrix using previously allocated storage. + * @param localMat view into the local matrix to populate + */ + virtual void extract( CRSMatrixView< real64, globalIndex > const & localMat ) const = 0; + + /** + * @brief Extract local part of the matrix using previously allocated storage and structure. + * @param localMat view into the local matrix to populate + */ + virtual void extract( CRSMatrixView< real64, globalIndex const > const & localMat ) const = 0; + + /** + * @brief Extract local part of the matrix + * @return the local matrix + */ + CRSMatrix< real64, globalIndex > extract() const + { + array1d< localIndex > rowLengths( numLocalRows() ); + getRowLengths( rowLengths.toView() ); + CRSMatrix< real64, globalIndex > localMat; + localMat.resizeFromRowCapacities< parallelHostPolicy >( numLocalRows(), numGlobalCols(), rowLengths.data() ); + extract( localMat.toView() ); + return localMat; + } + + /** + * @brief Extract block of the matrix corresponding to locally owned columns using previously allocated storage + * @param localMat the local matrix + * @note the column indices are shifted linearly to the local range + */ + virtual void extractLocal( CRSMatrixView< real64, localIndex > const & localMat ) const = 0; + + /** + * @brief Extract block of the matrix corresponding to locally owned columns + * @return the local matrix + * @note the column indices are shifted linearly to the local range + */ + CRSMatrix< real64, localIndex > extractLocal() const + { + array1d< localIndex > rowLengths( numLocalRows() ); + getRowLocalLengths( rowLengths.toView() ); + CRSMatrix< real64, localIndex > localMat; + localMat.resizeFromRowCapacities< parallelHostPolicy >( numLocalRows(), numGlobalCols(), rowLengths.data() ); + extractLocal( localMat.toView() ); + return localMat; + } + /** * @brief Populate a vector with row sums of @p this. * @param dst the target vector, must have the same row partitioning as @p this @@ -856,18 +938,6 @@ class MatrixBase : public virtual LinearOperator< VECTOR > */ virtual globalIndex jupper() const = 0; - /** - * @brief Returns the number of nonzeros in the local portion of the matrix - * @return the number of nonzeros in the local portion of the matrix - */ - virtual localIndex numLocalNonzeros() const = 0; - - /** - * @brief Returns the total number of nonzeros in the matrix - * @return the total number of nonzeros in the matrix - */ - virtual globalIndex numGlobalNonzeros() const = 0; - /** * @brief Returns the infinity norm of the matrix. * @return the value of infinity norm diff --git a/src/coreComponents/linearAlgebra/interfaces/direct/SuiteSparse.cpp b/src/coreComponents/linearAlgebra/interfaces/direct/SuiteSparse.cpp index a338ab2ae20..c7e525fb43a 100644 --- a/src/coreComponents/linearAlgebra/interfaces/direct/SuiteSparse.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/direct/SuiteSparse.cpp @@ -20,6 +20,7 @@ #include "codingUtilities/Utilities.hpp" #include "common/Stopwatch.hpp" +#include "common/TimingMacros.hpp" #include "linearAlgebra/common/common.hpp" #include "linearAlgebra/interfaces/InterfaceTypes.hpp" #include "linearAlgebra/utilities/Arnoldi.hpp" @@ -137,7 +138,7 @@ void factorize( SuiteSparseData & data, LinearSolverParameters const & params ) } // print the symbolic factorization - if( params.logLevel > 1 ) + if( params.logLevel >= 4 ) { umfpack_dl_report_symbolic( data.symbolic, data.control ); } @@ -159,7 +160,7 @@ void factorize( SuiteSparseData & data, LinearSolverParameters const & params ) } // print the numeric factorization - if( params.logLevel > 1 ) + if( params.logLevel >= 4 ) { umfpack_dl_report_numeric( data.symbolic, data.control ); } @@ -169,7 +170,7 @@ void setOptions( SuiteSparseData & data, LinearSolverParameters const & params ) { // Get the default control parameters umfpack_dl_defaults( data.control ); - data.control[UMFPACK_PRL] = params.logLevel > 1 ? 6 : 1; + data.control[UMFPACK_PRL] = params.logLevel; data.control[UMFPACK_ORDERING] = UMFPACK_ORDERING_BEST; } @@ -216,6 +217,7 @@ template< typename LAI > void SuiteSparse< LAI >::apply( Vector const & src, Vector & dst ) const { + GEOS_MARK_FUNCTION; doSolve( src, dst, false ); } @@ -271,7 +273,7 @@ void SuiteSparse< LAI >::solve( Vector const & rhs, condEst = estimateConditionNumberAdvanced(); if( m_result.residualReduction > condEst * precTol ) { - if( m_params.logLevel > 0 ) + if( m_params.logLevel > 0 && MpiWrapper::commRank( rhs.comm() ) ) { GEOS_WARNING( "SuiteSparse: failed to reduce residual below tolerance.\n" "Condition number estimate: " << condEst ); @@ -300,36 +302,45 @@ void SuiteSparse< LAI >::doSolve( Vector const & b, Vector & x, bool transpose ) GEOS_LAI_ASSERT_EQ( b.localSize(), x.localSize() ); GEOS_LAI_ASSERT_EQ( b.localSize(), matrix().numLocalRows() ); - m_export->exportVector( b, m_data->rhs ); + { + GEOS_MARK_SCOPE( export ); + m_export->exportVector( b, m_data->rhs ); + } - if( MpiWrapper::commRank( b.comm() ) == m_workingRank ) { - m_data->rhs.move( hostMemorySpace, false ); - m_data->sol.move( hostMemorySpace, true ); - - // To be able to use UMFPACK direct solver we need to disable floating point exceptions - LvArray::system::FloatingPointExceptionGuard guard; - - // Note: UMFPACK expects column-sparse matrix, but we have row-sparse, so we flip the transpose flag - SSlong const status = umfpack_dl_solve( transpose ? UMFPACK_A : UMFPACK_At, - m_data->rowPtr.data(), - m_data->colIndices.data(), - m_data->values.data(), - m_data->sol.data(), - m_data->rhs.data(), - m_data->numeric, - m_data->control, - m_data->info ); - - if( status < 0 ) + GEOS_MARK_SCOPE( solve ); + if( MpiWrapper::commRank( b.comm() ) == m_workingRank ) { - umfpack_dl_report_info( m_data->control, m_data->info ); - umfpack_dl_report_status( m_data->control, status ); - GEOS_ERROR( "SuiteSparse interface: umfpack_dl_solve failed." ); + m_data->rhs.move( hostMemorySpace, false ); + m_data->sol.move( hostMemorySpace, true ); + + // To be able to use UMFPACK direct solver we need to disable floating point exceptions + LvArray::system::FloatingPointExceptionGuard guard; + + // Note: UMFPACK expects column-sparse matrix, but we have row-sparse, so we flip the transpose flag + SSlong const status = umfpack_dl_solve( transpose ? UMFPACK_A : UMFPACK_At, + m_data->rowPtr.data(), + m_data->colIndices.data(), + m_data->values.data(), + m_data->sol.data(), + m_data->rhs.data(), + m_data->numeric, + m_data->control, + m_data->info ); + + if( status < 0 ) + { + umfpack_dl_report_info( m_data->control, m_data->info ); + umfpack_dl_report_status( m_data->control, status ); + GEOS_ERROR( "SuiteSparse interface: umfpack_dl_solve failed." ); + } } } - m_export->importVector( m_data->sol, x ); + { + GEOS_MARK_SCOPE( import ); + m_export->importVector( m_data->sol, x ); + } } template< typename LAI > @@ -357,7 +368,7 @@ real64 SuiteSparse< LAI >::estimateConditionNumberAdvanced() const GEOS_LAI_ASSERT( ready() ); localIndex constexpr numIterations = 4; - NormalOperator< LAI > const normalOperator( matrix() ); + NormalOperator< Matrix > const normalOperator( matrix() ); real64 const lambdaDirect = ArnoldiLargestEigenvalue( normalOperator, numIterations ); InverseNormalOperator< LAI, SuiteSparse > const inverseNormalOperator( matrix(), *this ); diff --git a/src/coreComponents/linearAlgebra/interfaces/direct/SuperLUDist.cpp b/src/coreComponents/linearAlgebra/interfaces/direct/SuperLUDist.cpp index 3e510657c59..1060029f4e4 100644 --- a/src/coreComponents/linearAlgebra/interfaces/direct/SuperLUDist.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/direct/SuperLUDist.cpp @@ -21,6 +21,7 @@ #include "codingUtilities/Utilities.hpp" #include "common/MpiWrapper.hpp" #include "common/Stopwatch.hpp" +#include "common/TimingMacros.hpp" #include "linearAlgebra/common/common.hpp" #include "linearAlgebra/interfaces/InterfaceTypes.hpp" #include "linearAlgebra/utilities/Arnoldi.hpp" @@ -112,7 +113,6 @@ struct SuperLUDistData array1d< int_t > rowPtr{}; ///< row pointers array1d< int_t > colIndices{}; ///< column indices array1d< double > values{}; ///< values - array1d< double > rhs{}; ///< rhs/solution vector values SuperMatrix mat{}; ///< SuperLU_Dist matrix format dScalePermstruct_t scalePerm{}; ///< data structure to scale and permute the matrix dLUstruct_t lu{}; ///< data structure to store the LU factorization @@ -129,7 +129,6 @@ struct SuperLUDistData rowPtr.resize( numLocalRows + 1 ); colIndices.resize( numLocalNonzeros ); values.resize( numLocalNonzeros ); - rhs.resize( numLocalRows ); dScalePermstructInit( numGlobalRows, numGlobalRows, &scalePerm ); dLUstructInit( numGlobalRows, &lu ); PStatInit( &stat ); @@ -211,6 +210,8 @@ template< typename LAI > void SuperLUDist< LAI >::apply( Vector const & src, Vector & dst ) const { + GEOS_MARK_FUNCTION; + GEOS_LAI_ASSERT( ready() ); GEOS_LAI_ASSERT( src.ready() ); GEOS_LAI_ASSERT( dst.ready() ); @@ -220,10 +221,10 @@ void SuperLUDist< LAI >::apply( Vector const & src, // To be able to use SuperLU_Dist solver we need to disable floating point exceptions LvArray::system::FloatingPointExceptionGuard guard; - // Export the rhs to a host-based array (this is required when vector is on GPU) - typename Matrix::Export vecExport; - vecExport.exportVector( src, m_data->rhs ); - m_data->rhs.move( hostMemorySpace, true ); + // pdgssvx operates on rhs/sol vector in-place + arrayView1d< real64 > const solution = dst.open(); + solution.setValues< serialPolicy >( src.values() ); + solution.move( hostMemorySpace, true ); // Call the linear equation solver to solve the matrix. real64 berr = 0.0; @@ -233,8 +234,8 @@ void SuperLUDist< LAI >::apply( Vector const & src, pdgssvx( &m_data->options, &m_data->mat, &m_data->scalePerm, - m_data->rhs.data(), - m_data->rhs.size(), + solution.data(), + solution.size(), 1, &m_data->grid, &m_data->lu, @@ -247,8 +248,7 @@ void SuperLUDist< LAI >::apply( Vector const & src, GEOS_LAI_ASSERT( !std::isnan( berr ) ); GEOS_LAI_ASSERT( !std::isinf( berr ) ); - // Import the solution back into the vector - vecExport.importVector( m_data->rhs, dst ); + dst.close(); } template< typename LAI > @@ -319,7 +319,7 @@ void SuperLUDist< LAI >::setOptions() { // Initialize options. set_default_options_dist( &m_data->options ); - m_data->options.PrintStat = m_params.logLevel > 1 ? YES : NO; + m_data->options.PrintStat = m_params.logLevel >= 3 ? YES : NO; m_data->options.Equil = m_params.direct.equilibrate ? YES : NO; m_data->options.ColPerm = getColPermType( m_params.direct.colPerm ); m_data->options.RowPerm = getRowPermType( m_params.direct.rowPerm ); @@ -327,7 +327,7 @@ void SuperLUDist< LAI >::setOptions() m_data->options.ReplaceTinyPivot = m_params.direct.replaceTinyPivot ? YES : NO; m_data->options.IterRefine = m_params.direct.iterativeRefine ? SLU_DOUBLE : NOREFINE; - if( m_params.logLevel > 0 ) + if( m_params.logLevel > 3 && MpiWrapper::commRank( m_data->grid.comm ) == 0 ) { print_sp_ienv_dist( &m_data->options ); print_options_dist( &m_data->options ); @@ -423,7 +423,7 @@ real64 SuperLUDist< LAI >::estimateConditionNumberAdvanced() const GEOS_LAI_ASSERT( ready() ); localIndex constexpr numIterations = 4; - NormalOperator< LAI > const normalOperator( matrix() ); + NormalOperator< Matrix > const normalOperator( matrix() ); real64 const lambdaDirect = ArnoldiLargestEigenvalue( normalOperator, numIterations ); InverseNormalOperator< LAI, SuperLUDist > const inverseNormalOperator( matrix(), *this ); diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreExport.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreExport.cpp index 7e335a47e1c..ef3b2444752 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreExport.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreExport.cpp @@ -157,10 +157,9 @@ void HypreExport::exportCRS( HypreMatrix const & mat, // Sort the values by column index after copying (some solvers expect this) forAll< hypre::execPolicy >( numRow, [rowOffsets, colIndices, values] GEOS_HYPRE_DEVICE ( HYPRE_Int const i ) { - using LvArray::sortedArrayManipulation::dualSort; - dualSort( colIndices.data() + rowOffsets[i], - colIndices.data() + rowOffsets[i + 1], - values.data() + rowOffsets[i] ); + LvArray::sortedArrayManipulation::dualSort( colIndices.data() + rowOffsets[i], + colIndices.data() + rowOffsets[i + 1], + values.data() + rowOffsets[i] ); } ); } @@ -175,7 +174,7 @@ void HypreExport::exportVector( HypreVector const & vec, // Gather vector on target rank, or just get the local part hypre_Vector * const localVector = m_targetRank < 0 ? hypre_ParVectorLocalVector( vec.unwrapped() ) - : (hypre_Vector *)hypre::parVectorToVectorAll( vec.unwrapped() ); + : (hypre_Vector *)hypre::parVectorToVector( vec.unwrapped(), m_targetRank ); GEOS_ERROR_IF( rank == m_targetRank && !localVector, "HypreExport: vector is empty on target rank" ); if( m_targetRank < 0 || m_targetRank == rank ) diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreKernels.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreKernels.hpp index eac1f2c468d..a85ff84ff9c 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreKernels.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreKernels.hpp @@ -209,72 +209,150 @@ makeSortedPermutation( HYPRE_Int const * const indices, } // namespace internal -template< typename SRC_COLMAP, typename DST_COLMAP > -void addEntriesRestricted( hypre_CSRMatrix const * const src_mat, - SRC_COLMAP const src_colmap, - hypre_CSRMatrix * const dst_mat, - DST_COLMAP const dst_colmap, - real64 const scale ) +template< typename KERNEL > +void addMatrixEntries( hypre_ParCSRMatrix const * const src, + hypre_ParCSRMatrix * const dst, + real64 const scale ) { - GEOS_LAI_ASSERT( src_mat != nullptr ); - GEOS_LAI_ASSERT( dst_mat != nullptr ); - - CSRData< true > src{ src_mat }; - CSRData< false > dst{ dst_mat }; - GEOS_LAI_ASSERT_EQ( src.nrow, dst.nrow ); - - if( src.ncol == 0 || isZero( scale ) ) + GEOS_LAI_ASSERT( src != nullptr ); + GEOS_LAI_ASSERT( dst != nullptr ); + KERNEL::launch( hypre_ParCSRMatrixDiag( src ), + hypre::ops::identity< HYPRE_Int >, + hypre_ParCSRMatrixDiag( dst ), + hypre::ops::identity< HYPRE_Int >, + scale ); + if( hypre_CSRMatrixNumCols( hypre_ParCSRMatrixOffd( dst ) ) > 0 ) { - return; + HYPRE_BigInt const * const src_colmap = hypre::getOffdColumnMap( src ); + HYPRE_BigInt const * const dst_colmap = hypre::getOffdColumnMap( dst ); + KERNEL::launch( hypre_ParCSRMatrixOffd( src ), + [src_colmap] GEOS_HYPRE_DEVICE ( auto i ){ return src_colmap[i]; }, + hypre_ParCSRMatrixOffd( dst ), + [dst_colmap] GEOS_HYPRE_DEVICE ( auto i ){ return dst_colmap[i]; }, + scale ); } +} - // Allocate contiguous memory to store sorted column permutations of each row - array1d< HYPRE_Int > const src_permutation_arr( hypre_CSRMatrixNumNonzeros( src_mat ) ); - array1d< HYPRE_Int > const dst_permutation_arr( hypre_CSRMatrixNumNonzeros( dst_mat ) ); - - arrayView1d< HYPRE_Int > const src_permutation = src_permutation_arr.toView(); - arrayView1d< HYPRE_Int > const dst_permutation = dst_permutation_arr.toView(); - // Each thread adds one row of src into dst - forAll< hypre::execPolicy >( dst.nrow, - [ src, - src_colmap, - dst, - dst_colmap, - scale, - src_permutation, - dst_permutation] GEOS_HYPRE_DEVICE ( HYPRE_Int const localRow ) +struct AddEntriesRestrictedKernel +{ + template< typename SRC_COLMAP, typename DST_COLMAP > + static void + launch( hypre_CSRMatrix const * const src_mat, + SRC_COLMAP const src_colmap, + hypre_CSRMatrix * const dst_mat, + DST_COLMAP const dst_colmap, + real64 const scale ) { - HYPRE_Int const src_offset = src.rowptr[localRow]; - HYPRE_Int const src_length = src.rowptr[localRow + 1] - src_offset; - HYPRE_Int const * const src_indices = src.colind + src_offset; - HYPRE_Real const * const src_values = src.values + src_offset; - HYPRE_Int * const src_perm = src_permutation.data() + src_offset; - - HYPRE_Int const dst_offset = dst.rowptr[localRow]; - HYPRE_Int const dst_length = dst.rowptr[localRow + 1] - dst_offset; - HYPRE_Int const * const dst_indices = dst.colind + dst_offset; - HYPRE_Real * const dst_values = dst.values + dst_offset; - HYPRE_Int * const dst_perm = dst_permutation.data() + dst_offset; - - // Since hypre does not store columns in sorted order, create a sorted "view" of src and dst rows - // TODO: it would be nice to cache the permutation arrays somewhere to avoid recomputing - internal::makeSortedPermutation( src_indices, src_length, src_perm, src_colmap ); - internal::makeSortedPermutation( dst_indices, dst_length, dst_perm, dst_colmap ); - - // Add entries looping through them in sorted column order, skipping src entries not in dst - for( HYPRE_Int i = 0, j = 0; i < dst_length && j < src_length; ++i ) + GEOS_LAI_ASSERT( src_mat != nullptr ); + GEOS_LAI_ASSERT( dst_mat != nullptr ); + + CSRData< true > src{ src_mat }; + CSRData< false > dst{ dst_mat }; + GEOS_LAI_ASSERT_EQ( src.nrow, dst.nrow ); + + if( src.ncol == 0 || isZero( scale ) ) + { + return; + } + + // Allocate contiguous memory to store sorted column permutations of each row + array1d< HYPRE_Int > const src_permutation_arr( src.nnz ); + array1d< HYPRE_Int > const dst_permutation_arr( dst.nnz ); + + arrayView1d< HYPRE_Int > const src_permutation = src_permutation_arr.toView(); + arrayView1d< HYPRE_Int > const dst_permutation = dst_permutation_arr.toView(); + // Each thread adds one row of src into dst + forAll< hypre::execPolicy >( dst.nrow, + [src, + src_colmap, + dst, + dst_colmap, + scale, + src_permutation, + dst_permutation ] GEOS_HYPRE_DEVICE ( HYPRE_Int const localRow ) { - while( j < src_length && src_colmap( src_indices[src_perm[j]] ) < dst_colmap( dst_indices[dst_perm[i]] ) ) + HYPRE_Int const src_offset = src.rowptr[localRow]; + HYPRE_Int const src_length = src.rowptr[localRow + 1] - src_offset; + HYPRE_Int const * const src_indices = src.colind + src_offset; + HYPRE_Real const * const src_values = src.values + src_offset; + HYPRE_Int * const src_perm = src_permutation.data() + src_offset; + + HYPRE_Int const dst_offset = dst.rowptr[localRow]; + HYPRE_Int const dst_length = dst.rowptr[localRow + 1] - dst_offset; + HYPRE_Int const * const dst_indices = dst.colind + dst_offset; + HYPRE_Real * const dst_values = dst.values + dst_offset; + HYPRE_Int * const dst_perm = dst_permutation.data() + dst_offset; + + // Since hypre does not store columns in sorted order, create a sorted "view" of src and dst rows + // TODO: it would be nice to cache the permutation arrays somewhere to avoid recomputing + internal::makeSortedPermutation( src_indices, src_length, src_perm, src_colmap ); + internal::makeSortedPermutation( dst_indices, dst_length, dst_perm, dst_colmap ); + + // Add entries looping through them in sorted column order, skipping src entries not in dst + for( HYPRE_Int i = 0, j = 0; i < dst_length && j < src_length; ++i ) { - ++j; + while( j < src_length && src_colmap( src_indices[src_perm[j]] ) < dst_colmap( dst_indices[dst_perm[i]] ) ) + { + ++j; + } + if( j < src_length && src_colmap( src_indices[src_perm[j]] ) == dst_colmap( dst_indices[dst_perm[i]] ) ) + { + dst_values[dst_perm[i]] += scale * src_values[src_perm[j++]]; + } } - if( j < src_length && src_colmap( src_indices[src_perm[j]] ) == dst_colmap( dst_indices[dst_perm[i]] ) ) + } ); + } +}; + +struct AddEntriesSamePatternKernel +{ + template< typename SRC_COLMAP, typename DST_COLMAP > + static void + launch( hypre_CSRMatrix const * const src_mat, + SRC_COLMAP const src_colmap, + hypre_CSRMatrix * const dst_mat, + DST_COLMAP const dst_colmap, + real64 const scale ) + { + GEOS_LAI_ASSERT( src_mat != nullptr ); + GEOS_LAI_ASSERT( dst_mat != nullptr ); + + CSRData< true > src{ src_mat }; + CSRData< false > dst{ dst_mat }; + GEOS_LAI_ASSERT_EQ( src.nrow, dst.nrow ); + + if( src.ncol == 0 || isZero( scale ) ) + { + return; + } + + // Each thread adds one row of src into dst + forAll< hypre::execPolicy >( dst.nrow, + [src, src_colmap, dst, dst_colmap, scale] GEOS_HYPRE_DEVICE ( HYPRE_Int const localRow ) + { + HYPRE_Int const src_offset = src.rowptr[localRow]; + HYPRE_Int const src_length = src.rowptr[localRow + 1] - src_offset; + HYPRE_Int const * const src_indices = src.colind + src_offset; + HYPRE_Real const * const src_values = src.values + src_offset; + + HYPRE_Int const dst_offset = dst.rowptr[localRow]; + HYPRE_Int const dst_length = dst.rowptr[localRow + 1] - dst_offset; + HYPRE_Int const * const dst_indices = dst.colind + dst_offset; + HYPRE_Real * const dst_values = dst.values + dst_offset; + + GEOS_ASSERT_EQ( src_offset, dst_offset ); + GEOS_ASSERT_EQ( src_length, dst_length ); + GEOS_DEBUG_VAR( src_length, dst_length, src_indices, dst_indices, src_colmap, dst_colmap ); + + // NOTE: this assumes that entries are in the exact same order, to avoid creating a sorted view + for( HYPRE_Int i = 0; i < dst_length; ++i ) { - dst_values[dst_perm[i]] += scale * src_values[src_perm[j++]]; + GEOS_ASSERT_EQ( src_colmap( src_indices[i] ), dst_colmap( dst_indices[i] ) ); + dst_values[i] += scale * src_values[i]; } - } - } ); -} + } ); + } +}; /// @endcond diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.cpp index c45c30a10b5..2abb59ffc2a 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.cpp @@ -689,17 +689,19 @@ void HypreMatrix::multiplyRAP( HypreMatrix const & R, dst.parCSRtoIJ( dst_parcsr ); } -void HypreMatrix::multiplyPtAP( HypreMatrix const & P, +void HypreMatrix::multiplyPtAP( HypreMatrix const & P1, + HypreMatrix const & P2, HypreMatrix & dst ) const { GEOS_LAI_ASSERT( ready() ); - GEOS_LAI_ASSERT( P.ready() ); - GEOS_LAI_ASSERT_EQ( numLocalRows(), P.numLocalRows() ); - GEOS_LAI_ASSERT_EQ( numLocalCols(), P.numLocalRows() ); + GEOS_LAI_ASSERT( P1.ready() ); + GEOS_LAI_ASSERT( P2.ready() ); + GEOS_LAI_ASSERT_EQ( numLocalRows(), P1.numLocalRows() ); + GEOS_LAI_ASSERT_EQ( numLocalCols(), P2.numLocalRows() ); - HYPRE_ParCSRMatrix const dst_parcsr = hypre_ParCSRMatrixRAPKT( P.unwrapped(), + HYPRE_ParCSRMatrix const dst_parcsr = hypre_ParCSRMatrixRAPKT( P1.unwrapped(), m_parcsr_mat, - P.unwrapped(), + P2.unwrapped(), 0 ); dst.parCSRtoIJ( dst_parcsr ); @@ -840,15 +842,12 @@ void HypreMatrix::separateComponentFilter( HypreMatrix & dst, integer const dofsPerNode ) const { GEOS_MARK_FUNCTION; + GEOS_LAI_ASSERT( ready() ); - localIndex const maxRowEntries = maxRowLength(); - integer const temp = maxRowEntries % dofsPerNode; - GEOS_LAI_ASSERT_EQ( temp, 0 ); - - CRSMatrix< real64 > tempMat; + CRSMatrix< real64, globalIndex > tempMat; - tempMat.resize( numLocalRows(), numGlobalCols(), maxRowEntries / dofsPerNode ); - CRSMatrixView< real64 > const tempMatView = tempMat.toView(); + tempMat.resize( numLocalRows(), numGlobalCols(), ( maxRowLengthLocal() + dofsPerNode - 1 ) / dofsPerNode ); + CRSMatrixView< real64, globalIndex > const tempMatView = tempMat.toView(); globalIndex const firstLocalRow = ilower(); globalIndex const firstLocalCol = jlower(); @@ -907,24 +906,14 @@ void HypreMatrix::addEntries( HypreMatrix const & src, { case MatrixPatternOp::Restrict: { - hypre::addEntriesRestricted( hypre_ParCSRMatrixDiag( src.unwrapped() ), - hypre::ops::identity< HYPRE_Int >, - hypre_ParCSRMatrixDiag( unwrapped() ), - hypre::ops::identity< HYPRE_Int >, - scale ); - if( hypre_CSRMatrixNumCols( hypre_ParCSRMatrixOffd( unwrapped() ) ) > 0 ) - { - HYPRE_BigInt const * const src_colmap = hypre::getOffdColumnMap( src.unwrapped() ); - HYPRE_BigInt const * const dst_colmap = hypre::getOffdColumnMap( unwrapped() ); - hypre::addEntriesRestricted( hypre_ParCSRMatrixOffd( src.unwrapped() ), - [src_colmap] GEOS_HYPRE_DEVICE ( auto i ) { return src_colmap[i]; }, - hypre_ParCSRMatrixOffd( unwrapped() ), - [dst_colmap] GEOS_HYPRE_DEVICE ( auto i ) { return dst_colmap[i]; }, - scale ); - } + hypre::addMatrixEntries< hypre::AddEntriesRestrictedKernel >( src.unwrapped(), unwrapped(), scale ); + break; + } + case MatrixPatternOp::Equal: + { + hypre::addMatrixEntries< hypre::AddEntriesSamePatternKernel >( src.unwrapped(), unwrapped(), scale ); break; } - case MatrixPatternOp::Same: case MatrixPatternOp::Subset: case MatrixPatternOp::Extend: { @@ -976,7 +965,7 @@ void HypreMatrix::clampEntries( real64 const lo, hypre::clampMatrixEntries( hypre_ParCSRMatrixOffd( m_parcsr_mat ), lo, hi, false ); } -localIndex HypreMatrix::maxRowLength() const +localIndex HypreMatrix::maxRowLengthLocal() const { GEOS_LAI_ASSERT( assembled() ); @@ -989,7 +978,7 @@ localIndex HypreMatrix::maxRowLength() const localMaxRowLength.max( (ia_diag[localRow + 1] - ia_diag[localRow]) + (ia_offd[localRow + 1] - ia_offd[localRow] ) ); } ); - return MpiWrapper::max( localMaxRowLength.get(), comm() ); + return localMaxRowLength.get(); } localIndex HypreMatrix::rowLength( globalIndex const globalRowIndex ) const @@ -1033,6 +1022,16 @@ void HypreMatrix::getRowLengths( arrayView1d< localIndex > const & lengths ) con } ); } +void HypreMatrix::getRowLocalLengths( arrayView1d< localIndex > const & lengths ) const +{ + GEOS_LAI_ASSERT( assembled() ); + HYPRE_Int const * const ia_diag = hypre_CSRMatrixI( hypre_ParCSRMatrixDiag( m_parcsr_mat ) ); + forAll< hypre::execPolicy >( numLocalRows(), [=] GEOS_HYPRE_DEVICE ( localIndex const localRow ) + { + lengths[localRow] = ia_diag[localRow + 1] - ia_diag[localRow]; + } ); +} + void HypreMatrix::getRowCopy( globalIndex const globalRowIndex, arraySlice1d< globalIndex > const & colIndices, arraySlice1d< real64 > const & values ) const @@ -1067,6 +1066,89 @@ void HypreMatrix::extractDiagonal( HypreVector & dst ) const dst.touch(); } +void HypreMatrix::extract( CRSMatrixView< real64, globalIndex > const & localMat ) const +{ + GEOS_LAI_ASSERT( ready() ); + GEOS_LAI_ASSERT_EQ( localMat.numRows(), numLocalRows() ); + GEOS_LAI_ASSERT_EQ( localMat.numColumns(), numGlobalCols() ); + + hypre::CSRData< true > const diag{ hypre_ParCSRMatrixDiag( unwrapped() ) }; + hypre::CSRData< true > const offd{ hypre_ParCSRMatrixOffd( unwrapped() ) }; + HYPRE_BigInt const * const colMap = hypre::getOffdColumnMap( unwrapped() ); + globalIndex const firstLocalCol = jlower(); + + forAll< hypre::execPolicy >( localMat.numRows(), + [localMat, diag, offd, + colMap, firstLocalCol] GEOS_HYPRE_DEVICE ( localIndex const localRow ) + { + localMat.removeNonZeros( localRow, localMat.getColumns( localRow ), localMat.numNonZeros( localRow ) ); + for( HYPRE_Int k = diag.rowptr[localRow]; k < diag.rowptr[localRow + 1]; ++k ) + { + globalIndex const col = firstLocalCol + diag.colind[k]; + localMat.insertNonZero( localRow, col, diag.values[k] ); + } + if( offd.ncol > 0 ) + { + for( HYPRE_Int k = offd.rowptr[localRow]; k < offd.rowptr[localRow + 1]; ++k ) + { + globalIndex const col = colMap[offd.colind[k]]; + localMat.insertNonZero( localRow, col, offd.values[k] ); + } + } + } ); +} + +void HypreMatrix::extract( CRSMatrixView< real64, globalIndex const > const & localMat ) const +{ + GEOS_LAI_ASSERT( ready() ); + GEOS_LAI_ASSERT_EQ( localMat.numRows(), numLocalRows() ); + GEOS_LAI_ASSERT_EQ( localMat.numColumns(), numGlobalCols() ); + + hypre::CSRData< true > const diag{ hypre_ParCSRMatrixDiag( unwrapped() ) }; + hypre::CSRData< true > const offd{ hypre_ParCSRMatrixOffd( unwrapped() ) }; + HYPRE_BigInt const * const colMap = hypre::getOffdColumnMap( unwrapped() ); + globalIndex const firstLocalCol = jlower(); + + localMat.zero(); + forAll< hypre::execPolicy >( localMat.numRows(), + [localMat, diag, offd, + colMap, firstLocalCol] GEOS_HYPRE_DEVICE ( localIndex const localRow ) + { + for( HYPRE_Int k = diag.rowptr[localRow]; k < diag.rowptr[localRow + 1]; ++k ) + { + globalIndex const col = firstLocalCol + diag.colind[k]; + localMat.addToRow< serialAtomic >( localRow, &col, &diag.values[k], 1 ); + } + if( offd.ncol > 0 ) + { + for( HYPRE_Int k = offd.rowptr[localRow]; k < offd.rowptr[localRow + 1]; ++k ) + { + globalIndex const col = colMap[offd.colind[k]]; + localMat.addToRow< serialAtomic >( localRow, &col, &offd.values[k], 1 ); + } + } + } ); +} + +void HypreMatrix::extractLocal( CRSMatrixView< real64, localIndex > const & localMat ) const +{ + GEOS_LAI_ASSERT( ready() ); + GEOS_LAI_ASSERT_EQ( localMat.numRows(), numLocalRows() ); + GEOS_LAI_ASSERT_EQ( localMat.numColumns(), numGlobalCols() ); + + hypre::CSRData< true > const diag{ hypre_ParCSRMatrixDiag( unwrapped() ) }; + forAll< hypre::execPolicy >( localMat.numRows(), + [localMat, diag] GEOS_HYPRE_DEVICE ( localIndex const localRow ) + { + localMat.removeNonZeros( localRow, localMat.getColumns( localRow ), localMat.numNonZeros( localRow ) ); + GEOS_ASSERT_GE( localMat.nonZeroCapacity( localRow ), diag.rowptr[localRow + 1] - diag.rowptr[localRow] ); + for( HYPRE_Int k = diag.rowptr[localRow]; k < diag.rowptr[localRow + 1]; ++k ) + { + localMat.insertNonZero( localRow, diag.colind[k], diag.values[k] ); + } + } ); +} + namespace { @@ -1315,17 +1397,21 @@ void HypreMatrix::write( string const & filename, { int const rank = MpiWrapper::commRank( comm() ); + globalIndex const numRows = numGlobalRows(); + globalIndex const numCols = numGlobalCols(); + globalIndex const numNonzeros = numGlobalNonzeros(); + // Write MatrixMarket header if( rank == 0 ) { std::ofstream os( filename ); GEOS_ERROR_IF( !os, GEOS_FMT( "Unable to open file for writing: {}", filename ) ); os << "%%MatrixMarket matrix coordinate real general\n"; - os << GEOS_FMT( "{} {} {}\n", numGlobalRows(), numGlobalCols(), numGlobalNonzeros() ); + os << GEOS_FMT( "{} {} {}\n", numRows, numCols, numNonzeros ); } // Write matrix values - if( numGlobalRows() > 0 && numGlobalCols() > 0 ) + if( numRows > 0 && numCols > 0 ) { // Copy distributed parcsr matrix in a local CSR matrix on every process with at least one row // Warning: works for a parcsr matrix that is smaller than 2^31-1 @@ -1341,13 +1427,15 @@ void HypreMatrix::write( string const & filename, std::ofstream os( filename, std::ios_base::app ); GEOS_ERROR_IF( !os, GEOS_FMT( "Unable to open file for writing on rank {}: {}", rank, filename ) ); char str[64]; + int const width = static_cast< int >( std::log10( std::max( csr.nrow, csr.ncol ) ) ) + 1; for( HYPRE_Int i = 0; i < csr.nrow; i++ ) { for( HYPRE_Int k = csr.rowptr[i]; k < csr.rowptr[i + 1]; k++ ) { // MatrixMarket row/col indices are 1-based - GEOS_FMT_TO( str, sizeof( str ), "{} {} {:>28.16e}\n", i + 1, csr.colind[k] + 1, csr.values[k] ); + GEOS_FMT_TO( str, sizeof( str ), "{1:>{0}} {2:>{0}} {3:>24.16e}\n", + width, i + 1, csr.colind[k] + 1, csr.values[k] ); os << str; } } diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.hpp index 41a7610c31c..8ae9b804669 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.hpp @@ -128,6 +128,9 @@ class HypreMatrix final : public virtual LinearOperator< HypreVector >, using MatrixBase::setDofManager; using MatrixBase::dofManager; using MatrixBase::create; + using MatrixBase::extract; + using MatrixBase::extractLocal; + using MatrixBase::multiplyPtAP; virtual void create( CRSMatrixView< real64 const, globalIndex const > const & localMatrix, localIndex const numLocalColumns, @@ -250,8 +253,9 @@ class HypreMatrix final : public virtual LinearOperator< HypreVector >, HypreMatrix const & P, HypreMatrix & dst ) const override; - virtual void multiplyPtAP( HypreMatrix const & P, - HypreMatrix & dst ) const override; + virtual void multiplyPtAP( Matrix const & P1, + Matrix const & P2, + Matrix & dst ) const override; virtual void gemv( real64 const alpha, HypreVector const & x, @@ -292,20 +296,28 @@ class HypreMatrix final : public virtual LinearOperator< HypreVector >, bool const excludeDiag ) override; /** - * @copydoc MatrixBase::maxRowLength + * @copydoc MatrixBase::maxRowLengthLocal */ - virtual localIndex maxRowLength() const override; + virtual localIndex maxRowLengthLocal() const override; virtual localIndex rowLength( globalIndex const globalRowIndex ) const override; virtual void getRowLengths( arrayView1d< localIndex > const & lengths ) const override; + virtual void getRowLocalLengths( arrayView1d< localIndex > const & lengths ) const override; + virtual void getRowCopy( globalIndex globalRowIndex, arraySlice1d< globalIndex > const & colIndices, arraySlice1d< real64 > const & values ) const override; virtual void extractDiagonal( HypreVector & dst ) const override; + virtual void extract( CRSMatrixView< real64, globalIndex > const & localMat ) const override; + + virtual void extract( CRSMatrixView< real64, globalIndex const > const & localMat ) const override; + + virtual void extractLocal( CRSMatrixView< real64, localIndex > const & localMat ) const override; + virtual void getRowSums( HypreVector & dst, RowSumType const rowSumType ) const override; diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.cpp index e4124a87e95..0c17e94793c 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.cpp @@ -85,12 +85,15 @@ void createAMG( LinearSolverParameters const & params, HYPRE_Int logLevel = (params.logLevel == 2 || params.logLevel >= 4) ? 1 : 0; GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetTol( precond.ptr, 0.0 ) ); - GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxIter( precond.ptr, 1 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxIter( precond.ptr, params.amg.numCycles ) ); GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetPrintLevel( precond.ptr, logLevel ) ); GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumFunctions( precond.ptr, params.dofsPerNode ) ); // Set maximum number of multigrid levels (default 25) - GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxLevels( precond.ptr, LvArray::integerConversion< HYPRE_Int >( params.amg.maxLevels ) ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxLevels( precond.ptr, params.amg.maxLevels ) ); + + // Set coarse grid max size (coarsening will stop once this limit is reached) + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxCoarseSize( precond.ptr, params.amg.maxCoarseSize ) ); // Set type of cycle (1: V-cycle (default); 2: W-cycle) GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCycleType( precond.ptr, hypre::getAMGCycleType( params.amg.cycleType ) ) ); @@ -135,11 +138,13 @@ void createAMG( LinearSolverParameters const & params, // Set smoother to be used (other options available, see hypre's documentation) // (default "gaussSeidel", i.e. local symmetric Gauss-Seidel) - if( params.amg.smootherType == LinearSolverParameters::AMG::SmootherType::ilu0 || + if( params.amg.smootherType == LinearSolverParameters::AMG::SmootherType::ilu || params.amg.smootherType == LinearSolverParameters::AMG::SmootherType::ilut ) { GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetSmoothType( precond.ptr, 5 ) ); - GEOS_LAI_CHECK_ERROR( HYPRE_ILUSetType( precond.ptr, hypre::getILUType( params.amg.smootherType ) ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetILUType( precond.ptr, hypre::getILUType( params.amg.smootherType ) ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetSmoothNumLevels( precond.ptr, LvArray::integerConversion< HYPRE_Int >( params.amg.maxLevels ) ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetSmoothNumSweeps( precond.ptr, LvArray::integerConversion< HYPRE_Int >( params.amg.numSweeps ) ) ); } else { @@ -247,6 +252,12 @@ void createILU( LinearSolverParameters const & params, GEOS_LAI_CHECK_ERROR( HYPRE_ILUSetDropThreshold( precond.ptr, params.ifact.threshold ) ); } + // Disable RCM reordering to avoid problems with mechanics + if( params.dofsPerNode > 1 ) + { + GEOS_LAI_CHECK_ERROR( HYPRE_ILUSetLocalReordering( precond.ptr, 0 ) ); + } + precond.setup = HYPRE_ILUSetup; precond.solve = HYPRE_ILUSolve; precond.destroy = HYPRE_ILUDestroy; @@ -261,6 +272,15 @@ void createRelaxation( LinearSolverParameters const & params, precond.destroy = hypre::relaxationDestroy; } +void createChebyshev( LinearSolverParameters const & params, + HyprePrecWrapper & precond ) +{ + GEOS_LAI_CHECK_ERROR( hypre::chebyshevCreate( precond.ptr, params.chebyshev.order, params.chebyshev.eigNumIter ) ); + precond.setup = hypre::chebyshevSetup; + precond.solve = hypre::chebyshevSolve; + precond.destroy = hypre::chebyshevDestroy; +} + } // namespace HyprePreconditioner::HyprePreconditioner( LinearSolverParameters params ) @@ -300,12 +320,16 @@ void HyprePreconditioner::create( DofManager const * const dofManager ) case LinearSolverParameters::PreconditionerType::bgs: case LinearSolverParameters::PreconditionerType::sgs: case LinearSolverParameters::PreconditionerType::l1jacobi: - case LinearSolverParameters::PreconditionerType::chebyshev: case LinearSolverParameters::PreconditionerType::l1sgs: { createRelaxation( m_params, *m_precond ); break; } + case LinearSolverParameters::PreconditionerType::chebyshev: + { + createChebyshev( m_params, *m_precond ); + break; + } case LinearSolverParameters::PreconditionerType::amg: { createAMG( m_params, *m_nullSpace, *m_precond ); @@ -317,7 +341,7 @@ void HyprePreconditioner::create( DofManager const * const dofManager ) hypre::mgr::createMGR( m_params, dofManager, *m_precond, *m_mgrData ); break; } - case LinearSolverParameters::PreconditionerType::iluk: + case LinearSolverParameters::PreconditionerType::ilu: case LinearSolverParameters::PreconditionerType::ilut: { createILU( m_params, *m_precond ); diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreUtils.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreUtils.cpp index cd20aa4a653..78c4f2b4895 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreUtils.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreUtils.cpp @@ -23,6 +23,8 @@ #include <_hypre_parcsr_mv.h> #include <_hypre_parcsr_ls.h> +#include + namespace geos { @@ -45,6 +47,68 @@ HYPRE_Vector parVectorToVectorAll( HYPRE_ParVector const vec ) } } +HYPRE_Vector parVectorToVector( HYPRE_ParVector const vec, int const targetRank ) +{ + // If input vector on device, clone to host so MPI can access its data + hypre_ParVector * const hostVec = hypre_ParVectorMemoryLocation( vec ) == HYPRE_MEMORY_HOST + ? (hypre_ParVector *)vec + : hypre_ParVectorCloneDeep_v2( vec, HYPRE_MEMORY_HOST ); + + MPI_Comm const comm = hypre_ParVectorComm( hostVec ); + int const rank = MpiWrapper::commRank( comm ); + int const numProcs = MpiWrapper::commSize( comm ); + int constexpr tag = 1337; + + HYPRE_Int const globalSize = LvArray::integerConversion< HYPRE_Int >( hypre_ParVectorGlobalSize( hostVec ) ); + HYPRE_Int const firstIndex = LvArray::integerConversion< HYPRE_Int >( hypre_ParVectorFirstIndex( hostVec ) ); + + hypre_Vector const * localVec = hypre_ParVectorLocalVector( hostVec ); + HYPRE_Int const localSize = hypre_VectorSize( localVec ); + HYPRE_Real const * const localData = hypre_VectorData( localVec ); + + array1d< HYPRE_Int > offsets; + if( rank == targetRank ) + { + offsets.resize( numProcs + 1 ); + offsets[numProcs] = globalSize; + } + MpiWrapper::gather( &firstIndex, 1, offsets.data(), 1, targetRank, comm ); + + hypre_Vector * newVec{}; + + if( rank == targetRank ) + { + newVec = hypre_SeqVectorCreate( globalSize ); + hypre_SeqVectorInitialize_v2( newVec, HYPRE_MEMORY_HOST ); + + std::vector< MPI_Request > requests( numProcs, MPI_REQUEST_NULL ); + for( int i = 0; i < numProcs; ++i ) + { + HYPRE_Real * const data = hypre_VectorData( newVec ) + offsets[i]; + if( i != rank ) + { + MpiWrapper::iRecv( data, offsets[i + 1] - offsets[i], i, tag, comm, &requests[i] ); + } + else + { + std::copy( localData, localData + localSize, data ); + } + } + MpiWrapper::waitAll( numProcs, requests.data(), MPI_STATUSES_IGNORE ); + } + else + { + MpiWrapper::send( localData, localSize, targetRank, tag, comm ); + } + + if( hostVec != vec ) + { + hypre_ParVectorDestroy( hostVec ); + } + + return (HYPRE_Vector)newVec; +} + HYPRE_Int dummySetup( HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, @@ -145,7 +209,7 @@ HYPRE_Int relaxationSetup( HYPRE_Solver solver, { hypre_TFree( data->norms, hypre::memoryLocation ); } - hypre_ParCSRComputeL1Norms( A, normType, nullptr, &data->norms ); + GEOS_LAI_CHECK_ERROR( hypre_ParCSRComputeL1Norms( A, normType, nullptr, &data->norms ) ); } return 0; } @@ -164,6 +228,7 @@ HYPRE_Int relaxationDestroy( HYPRE_Solver solver ) { // Refer to RelaxationData doxygen above for explanation of reinterpret_cast RelaxationData * const data = reinterpret_cast< RelaxationData * >( solver ); + GEOS_ASSERT( data != nullptr ); if( data->norms ) { hypre_TFree( data->norms, hypre::memoryLocation ); @@ -172,6 +237,111 @@ HYPRE_Int relaxationDestroy( HYPRE_Solver solver ) return 0; } +/** + * @brief Holds data needed by hypre's Chebyshev smoothing. + * + * See RelaxationData for explanation of the way this structure is used. + */ +struct ChebyshevData +{ + HYPRE_Int order = 1; ///< Chebyshev order + HYPRE_Int eigNumIter = 10; ///< Number of eigenvalue estimation iterations + + HYPRE_Real * diag{}; ///< Scaled diagonal values extracted during setup + HYPRE_Real * coefs{}; ///< Precomputed coefficients + + // Temporary vectors required by Chebyshev solve + HypreVector vtemp; + HypreVector ztemp; + HypreVector ptemp; + HypreVector rtemp; +}; + +HYPRE_Int chebyshevCreate( HYPRE_Solver & solver, + HYPRE_Int const order, + HYPRE_Int const eigNumIter ) +{ + ChebyshevData * const data = new ChebyshevData; + data->order = order; + data->eigNumIter = eigNumIter; + solver = reinterpret_cast< HYPRE_Solver >( data ); + return 0; +} + +HYPRE_Int chebyshevSetup( HYPRE_Solver solver, + HYPRE_ParCSRMatrix A, + HYPRE_ParVector b, + HYPRE_ParVector x ) +{ + GEOS_UNUSED_VAR( b, x ); + + // Refer to ChebyshevData doxygen above for explanation of reinterpret_cast + ChebyshevData * const data = reinterpret_cast< ChebyshevData * >( solver ); + data->vtemp.create( hypre_ParCSRMatrixNumRows( A ), hypre_ParCSRMatrixComm( A ) ); + data->ztemp.create( hypre_ParCSRMatrixNumRows( A ), hypre_ParCSRMatrixComm( A ) ); + data->ptemp.create( hypre_ParCSRMatrixNumRows( A ), hypre_ParCSRMatrixComm( A ) ); + data->rtemp.create( hypre_ParCSRMatrixNumRows( A ), hypre_ParCSRMatrixComm( A ) ); + + HYPRE_Real max_eig, min_eig; + if( data->eigNumIter > 0 ) + { + GEOS_LAI_CHECK_ERROR( hypre_ParCSRMaxEigEstimateCG( A, 1, data->eigNumIter, &max_eig, &min_eig ) ); + } + else + { + GEOS_LAI_CHECK_ERROR( hypre_ParCSRMaxEigEstimate( A, 1, &max_eig, &min_eig ) ); + } + + return hypre_ParCSRRelax_Cheby_Setup( A, + max_eig, + min_eig, + 0.3, // fraction + data->order, + 1, // scale + 0, // variant + &data->coefs, + &data->diag ); +} + +HYPRE_Int chebyshevSolve( HYPRE_Solver solver, + HYPRE_ParCSRMatrix A, + HYPRE_ParVector b, + HYPRE_ParVector x ) +{ + // Refer to ChebyshevData doxygen above for explanation of reinterpret_cast + ChebyshevData * const data = reinterpret_cast< ChebyshevData * >( solver ); + + return hypre_ParCSRRelax_Cheby_Solve( A, + b, + data->diag, + data->coefs, + data->order, + 1, // scale + 0, // variant + x, + data->vtemp.unwrapped(), + data->ztemp.unwrapped(), + data->ptemp.unwrapped(), + data->rtemp.unwrapped() ); +} + +HYPRE_Int chebyshevDestroy( HYPRE_Solver solver ) +{ + // Refer to ChebyshevData doxygen above for explanation of reinterpret_cast + ChebyshevData * const data = reinterpret_cast< ChebyshevData * >( solver ); + GEOS_ASSERT( data != nullptr ); + if( data->diag ) + { + hypre_TFree( data->diag, hypre::memoryLocation ); + } + if( data->coefs ) + { + hypre_TFree( data->coefs, HYPRE_MEMORY_HOST ); + } + delete data; + return 0; +} + } // namespace hypre } // namespace geos diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreUtils.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreUtils.hpp index 8f162a85fea..ff72a2af3dc 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreUtils.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreUtils.hpp @@ -201,15 +201,25 @@ inline HYPRE_BigInt const * toHypreBigInt( geos::globalIndex const * const index } /** - * @brief Gather a parallel vector on a every rank. + * @brief Gather a parallel vector on every rank. * @param vec the vector to gather - * @return A newly allocated serial vector (may be null on ranks that don't have any elements) + * @return a newly allocated serial vector (may be null on ranks that don't have any elements) + * @note caller takes ownership and must dispose of the vector appropriately * * This is a wrapper around hypre_ParVectorToVectorAll() that works for both host-based * and device-based vectors without relying on Unified Memory. */ HYPRE_Vector parVectorToVectorAll( HYPRE_ParVector const vec ); +/** + * @brief Gather a parallel vector onto a single rank. + * @param vec the vector to gather + * @param targetRank the rank to gather the vector on + * @return a newly allocated serial vector on @p targetRank, nullptr on the rest + * @note caller takes ownership and must dispose of the vector appropriately + */ +HYPRE_Vector parVectorToVector( HYPRE_ParVector const vec, int const targetRank ); + /** * @brief Dummy function that does nothing but conform to hypre's signature for preconditioner setup/apply functions. * @return always 0 (success) @@ -283,6 +293,50 @@ HYPRE_Int relaxationSolve( HYPRE_Solver solver, */ HYPRE_Int relaxationDestroy( HYPRE_Solver solver ); +/** + * @brief Create a Chebyshev smoother. + * @param solver the solver + * @param order Chebyshev order (degree) + * @param eigNumIter number of eigenvalue estimation iterations + * @return always 0 + */ +HYPRE_Int chebyshevCreate( HYPRE_Solver & solver, + HYPRE_Int const order, + HYPRE_Int const eigNumIter ); + +/** + * @brief Setup a Chebyshev smoother. + * @param solver the solver + * @param A the matrix + * @param b the rhs vector (unused) + * @param x the solution vector (unused) + * @return hypre error code + */ +HYPRE_Int chebyshevSetup( HYPRE_Solver solver, + HYPRE_ParCSRMatrix A, + HYPRE_ParVector b, + HYPRE_ParVector x ); + +/** + * @brief Solve with a Chebyshev smoother. + * @param solver the solver + * @param A the matrix + * @param b the rhs vector (unused) + * @param x the solution vector (unused) + * @return hypre error code + */ +HYPRE_Int chebyshevSolve( HYPRE_Solver solver, + HYPRE_ParCSRMatrix A, + HYPRE_ParVector b, + HYPRE_ParVector x ); + +/** + * @brief Destroy a Chebyshev smoother. + * @param solver the solver + * @return always 0 + */ +HYPRE_Int chebyshevDestroy( HYPRE_Solver solver ); + /** * @brief Returns hypre's identifier of the AMG cycle type. * @param type AMG cycle type @@ -378,7 +432,7 @@ inline HYPRE_Int getILUType( LinearSolverParameters::AMG::SmootherType const typ { static map< LinearSolverParameters::AMG::SmootherType, HYPRE_Int > const typeMap = { - { LinearSolverParameters::AMG::SmootherType::ilu0, 0 }, + { LinearSolverParameters::AMG::SmootherType::ilu, 0 }, { LinearSolverParameters::AMG::SmootherType::ilut, 1 }, }; return findOption( typeMap, type, "ILU", "HyprePreconditioner" ); @@ -462,7 +516,7 @@ inline HYPRE_Int getILUType( LinearSolverParameters::PreconditionerType const ty { static map< LinearSolverParameters::PreconditionerType, HYPRE_Int > const typeMap = { - { LinearSolverParameters::PreconditionerType::iluk, 0 }, + { LinearSolverParameters::PreconditionerType::ilu, 0 }, { LinearSolverParameters::PreconditionerType::ilut, 1 }, }; return findOption( typeMap, type, "ILU", "HyprePreconditioner" ); diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreVector.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreVector.cpp index 7353fc698ce..76815ffaa7d 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreVector.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreVector.cpp @@ -379,7 +379,7 @@ void HypreVector::write( string const & filename, for( HYPRE_Int i = 0; i < size; i++ ) { - GEOS_FMT_TO( str, sizeof( str ), "{:>28.16e}\n", data[i] ); + GEOS_FMT_TO( str, sizeof( str ), "{:>24.16e}\n", data[i] ); os << str; } } diff --git a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.cpp b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.cpp index 3d2341cfab7..9541d0d9b0d 100644 --- a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.cpp @@ -659,13 +659,12 @@ void PetscMatrix::multiplyRAP( PetscMatrix const & R, GEOS_LAI_ASSERT( ready() ); GEOS_LAI_ASSERT( P.ready() ); GEOS_LAI_ASSERT( R.ready() ); - GEOS_LAI_ASSERT_EQ( R.numGlobalCols(), numGlobalRows() ); - GEOS_LAI_ASSERT_EQ( numGlobalCols(), P.numGlobalRows() ); + GEOS_LAI_ASSERT_EQ( R.numLocalCols(), numLocalRows() ); + GEOS_LAI_ASSERT_EQ( numLocalCols(), P.numLocalRows() ); dst.reset(); GEOS_LAI_CHECK_ERROR( MatMatMatMult( R.m_mat, m_mat, P.m_mat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &dst.m_mat ) ); - dst.m_assembled = true; } @@ -674,18 +673,14 @@ void PetscMatrix::multiplyPtAP( PetscMatrix const & P, { GEOS_LAI_ASSERT( ready() ); GEOS_LAI_ASSERT( P.ready() ); - GEOS_LAI_ASSERT_EQ( numGlobalRows(), P.numGlobalRows() ); - GEOS_LAI_ASSERT_EQ( numGlobalCols(), P.numGlobalRows() ); + GEOS_LAI_ASSERT_EQ( numLocalRows(), P.numLocalRows() ); + GEOS_LAI_ASSERT_EQ( numLocalCols(), P.numLocalRows() ); dst.reset(); // To be able to use the PtAP product in some cases, we need to disable floating point exceptions - { - LvArray::system::FloatingPointExceptionGuard guard( FE_ALL_EXCEPT ); - - GEOS_LAI_CHECK_ERROR( MatPtAP( m_mat, P.m_mat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &dst.m_mat ) ); - } - + LvArray::system::FloatingPointExceptionGuard guard( FE_ALL_EXCEPT ); + GEOS_LAI_CHECK_ERROR( MatPtAP( m_mat, P.m_mat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &dst.m_mat ) ); dst.m_assembled = true; } @@ -701,17 +696,16 @@ void PetscMatrix::transpose( PetscMatrix & dst ) const void PetscMatrix::separateComponentFilter( PetscMatrix & dst, integer const dofsPerNode ) const { - localIndex const maxRowEntries = maxRowLength(); - GEOS_LAI_ASSERT_EQ( maxRowEntries % dofsPerNode, 0 ); + GEOS_LAI_ASSERT( ready() ); - CRSMatrix< real64 > tempMat; - tempMat.resize( numLocalRows(), numGlobalCols(), maxRowEntries / dofsPerNode ); - CRSMatrixView< real64 > const tempMatView = tempMat.toView(); + CRSMatrix< real64, globalIndex > tempMat; + tempMat.resize( numLocalRows(), numGlobalCols(), ( maxRowLengthLocal() + dofsPerNode - 1 ) / dofsPerNode ); + CRSMatrixView< real64, globalIndex > const tempMatView = tempMat.toView(); PetscInt firstRow, lastRow; GEOS_LAI_CHECK_ERROR( MatGetOwnershipRange( m_mat, &firstRow, &lastRow ) ); - auto const getComponent = [dofsPerNode] ( auto i ) + auto const getComponent = [dofsPerNode] ( auto const i ) { return LvArray::integerConversion< integer >( i % dofsPerNode ); }; @@ -730,7 +724,7 @@ void PetscMatrix::separateComponentFilter( PetscMatrix & dst, tempMatView.insertNonZero( row - firstRow, cols[k], vals[k] ); } } - GEOS_LAI_CHECK_ERROR( MatRestoreRow( m_mat, row, &numEntries, nullptr, &vals ) ); + GEOS_LAI_CHECK_ERROR( MatRestoreRow( m_mat, row, &numEntries, &cols, &vals ) ); } dst.create( tempMatView.toViewConst(), numLocalCols(), MPI_COMM_GEOSX ); @@ -785,7 +779,7 @@ void PetscMatrix::addEntries( PetscMatrix const & src, switch( op ) { - case MatrixPatternOp::Same: + case MatrixPatternOp::Equal: { GEOS_LAI_CHECK_ERROR( MatAXPY( m_mat, scale, src.m_mat, SAME_NONZERO_PATTERN ) ); break; @@ -851,7 +845,7 @@ void PetscMatrix::clampEntries( real64 const lo, } ); } -localIndex PetscMatrix::maxRowLength() const +localIndex PetscMatrix::maxRowLengthLocal() const { GEOS_LAI_ASSERT( assembled() ); RAJA::ReduceMax< parallelHostReduce, localIndex > maxLocalLength( 0 ); @@ -860,16 +854,17 @@ localIndex PetscMatrix::maxRowLength() const { maxLocalLength.max( rowLength( globalRow ) ); } ); - return MpiWrapper::max( maxLocalLength.get(), comm() ); + return maxLocalLength.get(); } localIndex PetscMatrix::rowLength( globalIndex const globalRowIndex ) const { GEOS_LAI_ASSERT( assembled() ); PetscInt ncols; - GEOS_LAI_CHECK_ERROR( MatGetRow( m_mat, LvArray::integerConversion< PetscInt >( globalRowIndex ), &ncols, nullptr, nullptr ) ); - localIndex const nnz = ncols; - GEOS_LAI_CHECK_ERROR( MatRestoreRow( m_mat, LvArray::integerConversion< PetscInt >( globalRowIndex ), &ncols, nullptr, nullptr ) ); + PetscInt const row = LvArray::integerConversion< PetscInt >( globalRowIndex ); + GEOS_LAI_CHECK_ERROR( MatGetRow( m_mat, row, &ncols, nullptr, nullptr ) ); + localIndex const nnz = LvArray::integerConversion< localIndex >( ncols ); + GEOS_LAI_CHECK_ERROR( MatRestoreRow( m_mat, row, &ncols, nullptr, nullptr ) ); return nnz; } @@ -884,6 +879,29 @@ void PetscMatrix::getRowLengths( arrayView1d< localIndex > const & lengths ) con } ); } +void PetscMatrix::getRowLocalLengths( arrayView1d< localIndex > const & lengths ) const +{ + GEOS_LAI_ASSERT( assembled() ); + globalIndex const rowOffset = ilower(); + PetscInt firtsCol, lastCol; + GEOS_LAI_CHECK_ERROR( MatGetOwnershipRangeColumn( m_mat, &firtsCol, &lastCol ) ); + auto const isLocalColumn = [firtsCol, lastCol]( PetscInt const c ) + { + return firtsCol <= c && c < lastCol; + }; + // Can't use parallel policy here, because of PETSc's single row checkout policy + forAll< serialPolicy >( numLocalRows(), [=]( localIndex const localRow ) + { + PetscInt ncols; + PetscInt const * cols; + PetscInt const row = LvArray::integerConversion< PetscInt >( localRow + rowOffset ); + GEOS_LAI_CHECK_ERROR( MatGetRow( m_mat, row, &ncols, &cols, nullptr ) ); + auto const count = std::count_if( cols, cols + ncols, isLocalColumn ); + GEOS_LAI_CHECK_ERROR( MatRestoreRow( m_mat, row, &ncols, &cols, nullptr ) ); + lengths[localRow] = LvArray::integerConversion< localIndex >( count ); + } ); +} + void PetscMatrix::getRowCopy( globalIndex const globalRow, arraySlice1d< globalIndex > const & colIndices, arraySlice1d< real64 > const & values ) const @@ -892,17 +910,19 @@ void PetscMatrix::getRowCopy( globalIndex const globalRow, GEOS_LAI_ASSERT_GE( globalRow, ilower() ); GEOS_LAI_ASSERT_GT( iupper(), globalRow ); + PetscInt const row = LvArray::integerConversion< PetscInt >( globalRow ); + PetscScalar const * vals; PetscInt const * inds; PetscInt numEntries; - GEOS_LAI_CHECK_ERROR( MatGetRow( m_mat, LvArray::integerConversion< PetscInt >( globalRow ), &numEntries, &inds, &vals ) ); + GEOS_LAI_CHECK_ERROR( MatGetRow( m_mat, row, &numEntries, &inds, &vals ) ); GEOS_LAI_ASSERT_GE( colIndices.size(), numEntries ); GEOS_LAI_ASSERT_GE( values.size(), numEntries ); std::copy( inds, inds + numEntries, colIndices.dataIfContiguous() ); std::copy( vals, vals + numEntries, values.dataIfContiguous() ); - GEOS_LAI_CHECK_ERROR( MatRestoreRow( m_mat, LvArray::integerConversion< PetscInt >( globalRow ), &numEntries, &inds, &vals ) ); + GEOS_LAI_CHECK_ERROR( MatRestoreRow( m_mat, row, &numEntries, &inds, &vals ) ); } void PetscMatrix::extractDiagonal( PetscVector & dst ) const @@ -915,6 +935,98 @@ void PetscMatrix::extractDiagonal( PetscVector & dst ) const dst.touch(); } +void PetscMatrix::extract( CRSMatrixView< real64, globalIndex > const & localMat ) const +{ + GEOS_LAI_ASSERT( ready() ); + GEOS_LAI_ASSERT_EQ( localMat.numRows(), numLocalRows() ); + GEOS_LAI_ASSERT_EQ( localMat.numColumns(), numGlobalCols() ); + + PetscInt firstRow, lastRow; + GEOS_LAI_CHECK_ERROR( MatGetOwnershipRange( m_mat, &firstRow, &lastRow ) ); + + localMat.move( LvArray::MemorySpace::host, false ); + for( PetscInt row = firstRow; row < lastRow; ++row ) + { + PetscInt numEntries; + PetscInt const * cols; + PetscReal const * vals; + GEOS_LAI_CHECK_ERROR( MatGetRow( m_mat, row, &numEntries, &cols, &vals ) ); + localIndex const localRow = LvArray::integerConversion< localIndex >( row - firstRow ); + + localMat.removeNonZeros( localRow, localMat.getColumns( localRow ), localMat.numNonZeros( localRow ) ); + for( int k = 0; k < numEntries; ++k ) + { + localMat.insertNonZero( localRow, cols[k], vals[k] ); + } + GEOS_LAI_CHECK_ERROR( MatRestoreRow( m_mat, row, &numEntries, &cols, &vals ) ); + } +} + +void PetscMatrix::extract( CRSMatrixView< real64, globalIndex const > const & localMat ) const +{ + GEOS_LAI_ASSERT( ready() ); + GEOS_LAI_ASSERT_EQ( localMat.numRows(), numLocalRows() ); + GEOS_LAI_ASSERT_EQ( localMat.numColumns(), numGlobalCols() ); + + PetscInt firstRow, lastRow; + GEOS_LAI_CHECK_ERROR( MatGetOwnershipRange( m_mat, &firstRow, &lastRow ) ); + + localMat.zero(); + localMat.move( LvArray::MemorySpace::host, false ); + for( PetscInt row = firstRow; row < lastRow; ++row ) + { + PetscInt numEntries; + PetscInt const * cols; + PetscReal const * vals; + GEOS_LAI_CHECK_ERROR( MatGetRow( m_mat, row, &numEntries, &cols, &vals ) ); + localIndex const localRow = LvArray::integerConversion< localIndex >( row - firstRow ); + + for( int k = 0; k < numEntries; ++k ) + { + localMat.addToRow< serialAtomic >( localRow, &cols[k], &vals[k], 1 ); + } + GEOS_LAI_CHECK_ERROR( MatRestoreRow( m_mat, row, &numEntries, &cols, &vals ) ); + } +} + +void PetscMatrix::extractLocal( CRSMatrixView< real64, localIndex > const & localMat ) const +{ + GEOS_LAI_ASSERT( ready() ); + GEOS_LAI_ASSERT_EQ( localMat.numRows(), numLocalRows() ); + GEOS_LAI_ASSERT_EQ( localMat.numColumns(), numGlobalCols() ); + + PetscInt firstRow, lastRow; + GEOS_LAI_CHECK_ERROR( MatGetOwnershipRange( m_mat, &firstRow, &lastRow ) ); + + PetscInt firstCol, lastCol; + GEOS_LAI_CHECK_ERROR( MatGetOwnershipRangeColumn( m_mat, &firstCol, &lastCol ) ); + auto const isLocalColumn = [firstCol, lastCol]( PetscInt const c ) + { + return firstCol <= c && c < lastCol; + }; + + localMat.move( LvArray::MemorySpace::host, false ); + for( PetscInt row = firstRow; row < lastRow; ++row ) + { + PetscInt numEntries; + PetscInt const * cols; + PetscReal const * vals; + GEOS_LAI_CHECK_ERROR( MatGetRow( m_mat, row, &numEntries, &cols, &vals ) ); + localIndex const localRow = LvArray::integerConversion< localIndex >( row - firstRow ); + + localMat.removeNonZeros( localRow, localMat.getColumns( localRow ), localMat.numNonZeros( localRow ) ); + for( int k = 0; k < numEntries; ++k ) + { + if( isLocalColumn( cols[k] ) ) + { + localIndex const localCol = LvArray::integerConversion< localIndex >( cols[k] - firstCol ); + localMat.insertNonZero( localRow, localCol, vals[k] ); + } + } + GEOS_LAI_CHECK_ERROR( MatRestoreRow( m_mat, row, &numEntries, &cols, &vals ) ); + } +} + namespace { @@ -924,11 +1036,12 @@ double reduceRow( Mat mat, R reducer ) { PetscScalar const * vals; - PetscInt const * inds; + PetscInt const * cols; PetscInt numEntries; - GEOS_LAI_CHECK_ERROR( MatGetRow( mat, globalRow, &numEntries, &inds, &vals ) ); + PetscInt const row = LvArray::integerConversion< PetscInt >( globalRow ); + GEOS_LAI_CHECK_ERROR( MatGetRow( mat, row, &numEntries, &cols, &vals ) ); PetscScalar const res = std::accumulate( vals, vals + numEntries, 0.0, reducer ); - GEOS_LAI_CHECK_ERROR( MatRestoreRow( mat, LvArray::integerConversion< PetscInt >( globalRow ), &numEntries, &inds, &vals ) ); + GEOS_LAI_CHECK_ERROR( MatRestoreRow( mat, row, &numEntries, &cols, &vals ) ); return res; } @@ -1287,4 +1400,4 @@ void PetscMatrix::write( string const & filename, GEOS_LAI_CHECK_ERROR( PetscViewerDestroy( &viewer ) ); } -} // end geosx namespace +} // end geos namespace diff --git a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.hpp b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.hpp index 1ebad306c2c..4010aecd9f1 100644 --- a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.hpp @@ -118,6 +118,9 @@ class PetscMatrix final : public virtual LinearOperator< PetscVector >, using MatrixBase::residual; using MatrixBase::setDofManager; using MatrixBase::dofManager; + using MatrixBase::extract; + using MatrixBase::extractLocal; + using MatrixBase::multiplyPtAP; virtual void createWithLocalSize( localIndex const localRows, localIndex const localCols, @@ -280,14 +283,22 @@ class PetscMatrix final : public virtual LinearOperator< PetscVector >, /** * @copydoc MatrixBase::maxRowLength */ - virtual localIndex maxRowLength() const override; + virtual localIndex maxRowLengthLocal() const override; virtual localIndex rowLength( globalIndex const globalRowIndex ) const override; virtual void getRowLengths( arrayView1d< localIndex > const & lengths ) const override; + virtual void getRowLocalLengths( arrayView1d< localIndex > const & lengths ) const override; + virtual void extractDiagonal( PetscVector & dst ) const override; + virtual void extract( CRSMatrixView< real64, globalIndex > const & localMat ) const override; + + virtual void extract( CRSMatrixView< real64, globalIndex const > const & localMat ) const override; + + virtual void extractLocal( CRSMatrixView< real64, localIndex > const & localMat ) const override; + virtual void getRowSums( PetscVector & dst, RowSumType const rowSumType ) const override; diff --git a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscPreconditioner.cpp b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscPreconditioner.cpp index 3ea2a64e951..c98e5a40c9c 100644 --- a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscPreconditioner.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscPreconditioner.cpp @@ -60,7 +60,7 @@ PCType getPetscSmootherType( LinearSolverParameters::PreconditionerType const & { static std::map< LinearSolverParameters::PreconditionerType, PCType > const typeMap = { - { LinearSolverParameters::PreconditionerType::iluk, PCILU }, + { LinearSolverParameters::PreconditionerType::ilu, PCILU }, { LinearSolverParameters::PreconditionerType::ic, PCICC }, { LinearSolverParameters::PreconditionerType::jacobi, PCJACOBI }, { LinearSolverParameters::PreconditionerType::l1jacobi, PCJACOBI }, @@ -166,6 +166,9 @@ void createPetscAMG( LinearSolverParameters const & params, // Set max number of levels GEOS_LAI_CHECK_ERROR( PCGAMGSetNlevels( precond, params.amg.maxLevels ) ); + // Set coarse grid max size (coarsening will stop once this limit is reached) + GEOS_LAI_CHECK_ERROR( PCGAMGSetCoarseEqLim( precond, params.amg.maxCoarseSize ) ); + // TODO: need someone familiar with PETSc to take a look at this #if 0 GEOS_LAI_CHECK_ERROR( PCSetType( precond, PCHMG ) ); @@ -325,7 +328,7 @@ void PetscPreconditioner::setup( PetscMatrix const & mat ) case LinearSolverParameters::PreconditionerType::fgs: case LinearSolverParameters::PreconditionerType::bgs: case LinearSolverParameters::PreconditionerType::sgs: - case LinearSolverParameters::PreconditionerType::iluk: + case LinearSolverParameters::PreconditionerType::ilu: case LinearSolverParameters::PreconditionerType::ic: { createPetscSmoother( m_params, m_precond ); diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.cpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.cpp index 46dc04a9cf1..222983ad98a 100644 --- a/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.cpp @@ -29,6 +29,7 @@ #include #include #include +#include #include @@ -453,8 +454,8 @@ void EpetraMatrix::multiplyRAP( EpetraMatrix const & R, GEOS_LAI_ASSERT( ready() ); GEOS_LAI_ASSERT( R.ready() ); GEOS_LAI_ASSERT( P.ready() ); - GEOS_LAI_ASSERT_EQ( numGlobalRows(), R.numGlobalCols() ); - GEOS_LAI_ASSERT_EQ( numGlobalCols(), P.numGlobalRows() ); + GEOS_LAI_ASSERT_EQ( R.numLocalCols(), numLocalRows() ); + GEOS_LAI_ASSERT_EQ( numLocalCols(), P.numLocalRows() ); Epetra_CrsMatrix * result = nullptr; GEOS_LAI_CHECK_ERROR( ML_Epetra::ML_Epetra_RAP( unwrapped(), P.unwrapped(), R.unwrapped(), result, false ) ); @@ -469,21 +470,22 @@ void EpetraMatrix::multiplyRAP( EpetraMatrix const & R, void EpetraMatrix::multiplyPtAP( EpetraMatrix const & P, EpetraMatrix & dst ) const { - // TODO: ML_Epetra_PtAP does not work with long long indices, find a workaround? -#if 0 GEOS_LAI_ASSERT( ready() ); GEOS_LAI_ASSERT( P.ready() ); GEOS_LAI_ASSERT_EQ( numGlobalRows(), P.numGlobalRows() ); GEOS_LAI_ASSERT_EQ( numGlobalCols(), P.numGlobalRows() ); Epetra_CrsMatrix * result = nullptr; - GEOS_LAI_CHECK_ERROR( ML_Epetra::ML_Epetra_PtAP( unwrapped(), P.unwrapped(), result, false ) ); + // TODO: ML_Epetra_PtAP does not work with long long indices, find a workaround? +#if 0 + GEOS_LAI_CHECK_ERROR( ML_Epetra::ML_Epetra_PtAP( unwrapped(), P.unwrapped(), result ) ); +#else + GEOS_LAI_CHECK_ERROR( ML_Epetra::Epetra_PtAP( unwrapped(), P.unwrapped(), result ) ); +#endif + // After we switch to Epetra_CrsMatrix for storage, can avoid this copy dst.create( *result ); delete result; -#else - MatrixBase::multiplyPtAP( P, dst ); -#endif } void EpetraMatrix::gemv( real64 const alpha, @@ -543,20 +545,19 @@ void EpetraMatrix::transpose( EpetraMatrix & dst ) const void EpetraMatrix::separateComponentFilter( EpetraMatrix & dst, integer const dofsPerNode ) const { - localIndex const maxRowEntries = maxRowLength(); - GEOS_LAI_ASSERT_EQ( maxRowEntries % dofsPerNode, 0 ); + GEOS_LAI_ASSERT( ready() ); - CRSMatrix< real64 > tempMat; - tempMat.resize( numLocalRows(), numGlobalCols(), maxRowEntries / dofsPerNode ); - CRSMatrixView< real64 > const tempMatView = tempMat.toView(); + CRSMatrix< real64, globalIndex > tempMat; + tempMat.resize( numLocalRows(), numGlobalCols(), ( maxRowLengthLocal() + dofsPerNode - 1 ) / dofsPerNode ); globalIndex const firstLocalRow = ilower(); - auto const getComponent = [dofsPerNode] ( auto i ) + auto const getComponent = [dofsPerNode] ( auto const i ) { return LvArray::integerConversion< integer >( i % dofsPerNode ); }; - forAll< parallelHostPolicy >( numLocalRows(), [&] ( localIndex const localRow ) + forAll< parallelHostPolicy >( numLocalRows(), [this, getComponent, firstLocalRow, + tempMatView = tempMat.toView()] ( localIndex const localRow ) { int numEntries; int * columns; @@ -574,7 +575,7 @@ void EpetraMatrix::separateComponentFilter( EpetraMatrix & dst, } } ); - dst.create( tempMatView.toViewConst(), numLocalCols(), MPI_COMM_GEOSX ); + dst.create( tempMat.toViewConst(), numLocalCols(), MPI_COMM_GEOSX ); dst.setDofManager( dofManager() ); } @@ -615,39 +616,91 @@ real64 EpetraMatrix::clearRow( globalIndex const globalRow, namespace { -void addEntriesRestricted( Epetra_CrsMatrix const & src, - Epetra_CrsMatrix const & dst, +template< typename MAP > +void makeSortedPermutation( int const * const indices, + int const size, + int * const perm, + MAP map ) +{ + for( int i = 0; i < size; ++i ) + { + perm[i] = i; // std::iota + } + auto const comp = [indices, map] ( int i, int j ){ return map( indices[i] ) < map( indices[j] ); }; + LvArray::sortedArrayManipulation::makeSorted( perm, perm + size, comp ); +} + +struct CSRData +{ + int * rowptr{}; + int * colind{}; + double * values{}; + int nrow; + int ncol; + int nnz; + + explicit CSRData( Epetra_CrsMatrix const & mat ) + : nrow( mat.NumMyRows() ), + ncol( mat.NumMyCols() ), + nnz( mat.NumMyNonzeros() ) + { + mat.ExtractCrsDataPointers( rowptr, colind, values ); + } +}; + +void addEntriesRestricted( Epetra_CrsMatrix const & src_mat, + Epetra_CrsMatrix const & dst_mat, real64 const scale ) { - GEOS_LAI_ASSERT( src.NumMyRows() == dst.NumMyRows() ); + GEOS_LAI_ASSERT( src_mat.NumMyRows() == dst_mat.NumMyRows() ); + + CSRData src{ src_mat }; + CSRData dst{ dst_mat }; - if( isZero( scale ) ) + if( src.ncol == 0 || isZero( scale ) ) { return; } - int dst_length; - int * dst_indices; - double * dst_values; + array1d< int > const src_permutation( src.nnz ); + array1d< int > const dst_permutation( dst.nnz ); - int src_length; - int * src_indices; - double * src_values; - - for( int localRow = 0; localRow < dst.NumMyRows(); ++localRow ) + forAll< parallelHostPolicy >( dst.nrow, + [&src_mat, &dst_mat, src, dst, scale, + src_permutation = src_permutation.toView(), + dst_permutation = dst_permutation.toView()]( int const localRow ) { - dst.ExtractMyRowView( LvArray::integerConversion< int >( localRow ), dst_length, dst_values, dst_indices ); - src.ExtractMyRowView( LvArray::integerConversion< int >( localRow ), src_length, src_values, src_indices ); + + int const src_offset = src.rowptr[localRow]; + int const src_length = src.rowptr[localRow + 1] - src_offset; + int const * const src_indices = src.colind + src_offset; + double const * const src_values = src.values + src_offset; + int * const src_perm = src_permutation.data() + src_offset; + auto const src_colmap = [&src_mat]( int const i ) { return src_mat.GCID64( i ); }; + + int const dst_offset = dst.rowptr[localRow]; + int const dst_length = dst.rowptr[localRow + 1] - dst_offset; + int const * const dst_indices = dst.colind + dst_offset; + double * const dst_values = dst.values + dst_offset; + int * const dst_perm = dst_permutation.data() + dst_offset; + auto const dst_colmap = [&dst_mat]( int const i ) { return dst_mat.GCID64( i ); }; + + // Create a view of both matrix rows that is "sorted" w.r.t. GCID + makeSortedPermutation( src_indices, src_length, src_perm, src_colmap ); + makeSortedPermutation( dst_indices, dst_length, dst_perm, dst_colmap ); + for( int i = 0, j = 0; i < dst_length && j < src_length; ++i ) { - while( j < src_length && src_indices[j] < dst_indices[i] ) + while( j < src_length && src_colmap( src_indices[src_perm[j]] ) < dst_colmap( dst_indices[dst_perm[i]] ) ) + { ++j; - if( j < src_length && src_indices[j] == dst_indices[i] ) + } + if( j < src_length && src_colmap( src_indices[src_perm[j]] ) == dst_colmap( dst_indices[dst_perm[i]] ) ) { - dst_values[i] += scale * src_values[j++]; + dst_values[dst_perm[i]] += scale * src_values[src_perm[j++]]; } } - } + } ); } } // namespace @@ -663,7 +716,7 @@ void EpetraMatrix::addEntries( EpetraMatrix const & src, switch( op ) { - case MatrixPatternOp::Same: + case MatrixPatternOp::Equal: case MatrixPatternOp::Subset: { GEOS_LAI_CHECK_ERROR( EpetraExt::MatrixMatrix::Add( src.unwrapped(), false, scale, *m_matrix, 1.0 ) ); @@ -725,6 +778,12 @@ void EpetraMatrix::clampEntries( real64 const lo, } ); } +localIndex EpetraMatrix::maxRowLengthLocal() const +{ + GEOS_LAI_ASSERT( assembled() ); + return m_matrix->MaxNumEntries(); +} + localIndex EpetraMatrix::maxRowLength() const { GEOS_LAI_ASSERT( assembled() ); @@ -740,9 +799,30 @@ localIndex EpetraMatrix::rowLength( globalIndex const globalRowIndex ) const void EpetraMatrix::getRowLengths( arrayView1d< localIndex > const & lengths ) const { GEOS_LAI_ASSERT( assembled() ); - forAll< parallelHostPolicy >( numLocalRows(), [=]( localIndex const localRow ) + forAll< parallelHostPolicy >( numLocalRows(), [&mat = *m_matrix, lengths]( localIndex const localRow ) { - lengths[localRow] = m_matrix->NumMyEntries( LvArray::integerConversion< int >( localRow ) ); + lengths[localRow] = mat.NumMyEntries( LvArray::integerConversion< int >( localRow ) ); + } ); +} + +void EpetraMatrix::getRowLocalLengths( arrayView1d< localIndex > const & lengths ) const +{ + GEOS_LAI_ASSERT( assembled() ); + globalIndex const firstCol = jlower(); + globalIndex const lastCol = jupper(); + auto const isLocalColumn = [&mat = *m_matrix, firstCol, lastCol]( int const c ) + { + long long const gid = mat.GCID64( c ); + return firstCol <= gid && gid < lastCol; + }; + forAll< parallelHostPolicy >( numLocalRows(), [&mat = *m_matrix, lengths, isLocalColumn]( localIndex const localRow ) + { + int numEntries; + int * indicesPtr; + double * values; + GEOS_LAI_CHECK_ERROR( mat.ExtractMyRowView( localRow, numEntries, values, indicesPtr ) ); + auto const count = std::count_if( indicesPtr, indicesPtr + numEntries, isLocalColumn ); + lengths[localRow] = LvArray::integerConversion< localIndex >( count ); } ); } @@ -764,7 +844,7 @@ void EpetraMatrix::getRowCopy( globalIndex globalRow, GEOS_LAI_ASSERT_GE( values.size(), numEntries ); std::transform( indicesPtr, indicesPtr + numEntries, colIndices.begin(), - [&mat=*m_matrix]( int const c ){ return LvArray::integerConversion< globalIndex >( mat.GCID64( c ) ); } ); + [this]( int const c ){ return LvArray::integerConversion< globalIndex >( m_matrix->GCID64( c ) ); } ); std::copy( valuesPtr, valuesPtr + numEntries, values.begin() ); } @@ -778,6 +858,79 @@ void EpetraMatrix::extractDiagonal( EpetraVector & dst ) const dst.touch(); } +void EpetraMatrix::extract( CRSMatrixView< real64, globalIndex > const & localMat ) const +{ + GEOS_LAI_ASSERT( ready() ); + GEOS_LAI_ASSERT_EQ( localMat.numRows(), numLocalRows() ); + GEOS_LAI_ASSERT_EQ( localMat.numColumns(), numGlobalCols() ); + + forAll< parallelHostPolicy >( localMat.numRows(), [this, localMat] ( localIndex const localRow ) + { + int numEntries; + int * indicesPtr; + double * valuesPtr; + GEOS_LAI_CHECK_ERROR( m_matrix->ExtractMyRowView( localRow, numEntries, valuesPtr, indicesPtr ) ); + + localMat.removeNonZeros( localRow, localMat.getColumns( localRow ), localMat.numNonZeros( localRow ) ); + for( int k = 0; k < numEntries; ++k ) + { + globalIndex const col = m_matrix->GCID64( indicesPtr[k] ); + localMat.insertNonZero( localRow, col, valuesPtr[k] ); + } + } ); +} + +void EpetraMatrix::extract( CRSMatrixView< real64, globalIndex const > const & localMat ) const +{ + GEOS_LAI_ASSERT( ready() ); + GEOS_LAI_ASSERT_EQ( localMat.numRows(), numLocalRows() ); + GEOS_LAI_ASSERT_EQ( localMat.numColumns(), numGlobalCols() ); + + localMat.zero(); + forAll< parallelHostPolicy >( localMat.numRows(), [this, localMat] ( localIndex const localRow ) + { + int numEntries; + int * indicesPtr; + double * valuesPtr; + GEOS_LAI_CHECK_ERROR( m_matrix->ExtractMyRowView( localRow, numEntries, valuesPtr, indicesPtr ) ); + + for( int k = 0; k < numEntries; ++k ) + { + globalIndex const col = m_matrix->GCID64( indicesPtr[k] ); + localMat.addToRow< serialAtomic >( localRow, &col, &valuesPtr[k], 1 ); + } + } ); +} + +void EpetraMatrix::extractLocal( CRSMatrixView< real64, localIndex > const & localMat ) const +{ + GEOS_LAI_ASSERT( ready() ); + GEOS_LAI_ASSERT_EQ( localMat.numRows(), numLocalRows() ); + GEOS_LAI_ASSERT_EQ( localMat.numColumns(), numGlobalCols() ); + + globalIndex const firstCol = jlower(); + globalIndex const lastCol = jupper(); + + forAll< parallelHostPolicy >( localMat.numRows(), [=, & mat = *m_matrix] ( localIndex const localRow ) + { + int numEntries; + int * indicesPtr; + double * valuesPtr; + GEOS_LAI_CHECK_ERROR( mat.ExtractMyRowView( localRow, numEntries, valuesPtr, indicesPtr ) ); + + localMat.removeNonZeros( localRow, localMat.getColumns( localRow ), localMat.numNonZeros( localRow ) ); + for( int k = 0; k < numEntries; ++k ) + { + long long const globalCol = mat.GCID64( indicesPtr[k] ); + if( firstCol <= globalCol && globalCol < lastCol ) + { + localIndex const col = LvArray::integerConversion< localIndex >( globalCol - firstCol ); + localMat.insertNonZero( localRow, col, valuesPtr[k] ); + } + } + } ); +} + namespace { @@ -815,6 +968,7 @@ void rescaleRow( Epetra_CrsMatrix & mat, double * vals; mat.ExtractMyRowView( localRow, numEntries, vals ); double const scale = std::accumulate( vals, vals + numEntries, 0.0, reducer ); + GEOS_ASSERT_MSG( !isZero( scale ), "Zero row sum in row " << mat.GRID64( localRow ) ); std::transform( vals, vals + numEntries, vals, [scale]( double const v ){ return v / scale; } ); } diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.hpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.hpp index 44d907d37d5..ac9e9207104 100644 --- a/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.hpp @@ -107,6 +107,9 @@ class EpetraMatrix final : public virtual LinearOperator< EpetraVector >, using MatrixBase::residual; using MatrixBase::setDofManager; using MatrixBase::dofManager; + using MatrixBase::extract; + using MatrixBase::extractLocal; + using MatrixBase::multiplyPtAP; virtual void createWithLocalSize( localIndex const localRows, localIndex const localCols, @@ -266,6 +269,11 @@ class EpetraMatrix final : public virtual LinearOperator< EpetraVector >, real64 const hi, bool const excludeDiag ) override; + /** + * @copydoc MatrixBase::maxRowLengthLocal + */ + virtual localIndex maxRowLengthLocal() const override; + /** * @copydoc MatrixBase::maxRowLength */ @@ -275,12 +283,20 @@ class EpetraMatrix final : public virtual LinearOperator< EpetraVector >, virtual void getRowLengths( arrayView1d< localIndex > const & lengths ) const override; + virtual void getRowLocalLengths( arrayView1d< localIndex > const & lengths ) const override; + virtual void getRowCopy( globalIndex globalRow, arraySlice1d< globalIndex > const & colIndices, arraySlice1d< real64 > const & values ) const override; virtual void extractDiagonal( EpetraVector & dst ) const override; + virtual void extract( CRSMatrixView< real64, globalIndex > const & localMat ) const override; + + virtual void extract( CRSMatrixView< real64, globalIndex const > const & localMat ) const override; + + virtual void extractLocal( CRSMatrixView< real64, localIndex > const & localMat ) const override; + virtual void getRowSums( EpetraVector & dst, RowSumType const rowSumType ) const override; diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosPreconditioner.cpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosPreconditioner.cpp index 94963add498..3b4200c306e 100644 --- a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosPreconditioner.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosPreconditioner.cpp @@ -103,8 +103,8 @@ string getMLSmootherType( LinearSolverParameters::AMG::SmootherType const & valu { LinearSolverParameters::AMG::SmootherType::sgs, "symmetric Gauss-Seidel" }, { LinearSolverParameters::AMG::SmootherType::l1sgs, "symmetric Gauss-Seidel" }, { LinearSolverParameters::AMG::SmootherType::chebyshev, "Chebyshev" }, - { LinearSolverParameters::AMG::SmootherType::ic0, "IC" }, - { LinearSolverParameters::AMG::SmootherType::ilu0, "ILU" }, + { LinearSolverParameters::AMG::SmootherType::ic, "IC" }, + { LinearSolverParameters::AMG::SmootherType::ilu, "ILU" }, { LinearSolverParameters::AMG::SmootherType::ilut, "ILUT" }, }; @@ -154,14 +154,17 @@ createMLOperator( LinearSolverParameters const & params, list.set( "ML output", params.logLevel ); list.set( "max levels", params.amg.maxLevels ); + list.set( "coarse: max size", params.amg.maxCoarseSize ); list.set( "aggregation: type", "Uncoupled" ); list.set( "aggregation: threshold", params.amg.threshold ); - list.set( "PDE equations", params.dofsPerNode ); + list.set( "PDE equations", params.amg.numFunctions ); list.set( "smoother: sweeps", params.amg.numSweeps ); list.set( "prec type", getMLCycleType( params.amg.cycleType ) ); list.set( "smoother: type", getMLSmootherType( params.amg.smootherType ) ); list.set( "smoother: pre or post", getMLPreOrPostSmoothingType( params.amg.preOrPostSmoothing ) ); list.set( "coarse: type", getMLCoarseType( params.amg.coarseType ) ); + list.set( "repartition: enable", 1 ); + list.set( "repartition: partitioner", "ParMETIS" ); if( params.amg.nullSpaceType == LinearSolverParameters::AMG::NullSpaceType::constantModes ) { @@ -175,25 +178,24 @@ createMLOperator( LinearSolverParameters const & params, list.set( "null space: dimension", LvArray::integerConversion< integer >( nullSpacePointer.size( 0 ) ) ); } - std::unique_ptr< Epetra_Operator > precond = - std::make_unique< ML_Epetra::MultiLevelPreconditioner >( matrix, list ); - - return precond; + return std::make_unique< ML_Epetra::MultiLevelPreconditioner >( matrix, list, true ); } Ifpack::EPrecType getIfpackPrecondType( LinearSolverParameters::PreconditionerType const & type ) { static std::map< LinearSolverParameters::PreconditionerType, Ifpack::EPrecType > const typeMap = { - { LinearSolverParameters::PreconditionerType::iluk, Ifpack::ILU }, + { LinearSolverParameters::PreconditionerType::ilu, Ifpack::ILU }, { LinearSolverParameters::PreconditionerType::ilut, Ifpack::ILUT }, { LinearSolverParameters::PreconditionerType::ic, Ifpack::IC }, { LinearSolverParameters::PreconditionerType::ict, Ifpack::ICT }, { LinearSolverParameters::PreconditionerType::chebyshev, Ifpack::CHEBYSHEV }, { LinearSolverParameters::PreconditionerType::jacobi, Ifpack::POINT_RELAXATION }, + { LinearSolverParameters::PreconditionerType::l1jacobi, Ifpack::POINT_RELAXATION }, { LinearSolverParameters::PreconditionerType::fgs, Ifpack::POINT_RELAXATION }, { LinearSolverParameters::PreconditionerType::bgs, Ifpack::POINT_RELAXATION }, { LinearSolverParameters::PreconditionerType::sgs, Ifpack::POINT_RELAXATION }, + { LinearSolverParameters::PreconditionerType::l1sgs, Ifpack::POINT_RELAXATION }, { LinearSolverParameters::PreconditionerType::direct, Ifpack::AMESOS } }; @@ -202,11 +204,16 @@ Ifpack::EPrecType getIfpackPrecondType( LinearSolverParameters::PreconditionerTy } string getIfpackRelaxationType( LinearSolverParameters::PreconditionerType const & type ) { + // Trilinos does not implement l1 smoother variants, so we map them to regular + // smoothers for compatibility with input files designed for hypre. static std::map< LinearSolverParameters::PreconditionerType, string > const typeMap = { { LinearSolverParameters::PreconditionerType::jacobi, "Jacobi" }, + { LinearSolverParameters::PreconditionerType::l1jacobi, "Jacobi" }, { LinearSolverParameters::PreconditionerType::fgs, "Gauss-Seidel" }, + { LinearSolverParameters::PreconditionerType::bgs, "Gauss-Seidel" }, { LinearSolverParameters::PreconditionerType::sgs, "symmetric Gauss-Seidel" }, + { LinearSolverParameters::PreconditionerType::l1sgs, "symmetric Gauss-Seidel" }, }; GEOS_LAI_ASSERT_MSG( typeMap.count( type ) > 0, "Unsupported Trilinos/Ifpack preconditioner option: " << type ); @@ -245,6 +252,11 @@ createIfpackOperator( LinearSolverParameters const & params, list.set( "fact: ict level-of-fill", std::max( real64( params.ifact.fill ), 1.0 ) ); break; } + case Ifpack::CHEBYSHEV: + { + list.set( "relaxation: sweeps", params.chebyshev.order ); + list.set( "chebyshev: eigenvalue max iterations", params.chebyshev.eigNumIter ); + } default: { break; @@ -266,11 +278,19 @@ EpetraMatrix const & TrilinosPreconditioner::setupPreconditioningMatrix( EpetraM mat.separateComponentFilter( m_precondMatrix, m_params.dofsPerNode ); return m_precondMatrix; } - return mat; + else + { + // To avoid the issue of ML destructor crashing if matrix has been disposed of, + // we always perform setup on a copy of the input matrix + m_precondMatrix = mat; + } + return m_precondMatrix; } void TrilinosPreconditioner::setup( Matrix const & mat ) { + m_precond.reset(); + EpetraMatrix const & precondMat = setupPreconditioningMatrix( mat ); Base::setup( precondMat ); @@ -286,11 +306,13 @@ void TrilinosPreconditioner::setup( Matrix const & mat ) break; } case LinearSolverParameters::PreconditionerType::jacobi: + case LinearSolverParameters::PreconditionerType::l1jacobi: case LinearSolverParameters::PreconditionerType::fgs: case LinearSolverParameters::PreconditionerType::bgs: case LinearSolverParameters::PreconditionerType::sgs: + case LinearSolverParameters::PreconditionerType::l1sgs: case LinearSolverParameters::PreconditionerType::chebyshev: - case LinearSolverParameters::PreconditionerType::iluk: + case LinearSolverParameters::PreconditionerType::ilu: case LinearSolverParameters::PreconditionerType::ilut: case LinearSolverParameters::PreconditionerType::ic: case LinearSolverParameters::PreconditionerType::ict: @@ -301,7 +323,6 @@ void TrilinosPreconditioner::setup( Matrix const & mat ) } case LinearSolverParameters::PreconditionerType::none: { - m_precond.reset(); break; } default: diff --git a/src/coreComponents/linearAlgebra/solvers/BicgstabSolver.cpp b/src/coreComponents/linearAlgebra/solvers/BicgstabSolver.cpp index 8ba335f54dd..489567d2092 100644 --- a/src/coreComponents/linearAlgebra/solvers/BicgstabSolver.cpp +++ b/src/coreComponents/linearAlgebra/solvers/BicgstabSolver.cpp @@ -19,6 +19,7 @@ #include "BicgstabSolver.hpp" +#include "common/TimingMacros.hpp" #include "common/Stopwatch.hpp" #include "linearAlgebra/interfaces/InterfaceTypes.hpp" #include "common/LinearOperator.hpp" @@ -39,6 +40,7 @@ template< typename VECTOR > void BicgstabSolver< VECTOR >::solve( Vector const & b, Vector & x ) const { + GEOS_MARK_FUNCTION; Stopwatch watch; // Define vectors @@ -91,8 +93,8 @@ void BicgstabSolver< VECTOR >::solve( Vector const & b, // Compute r0.rk real64 const rho = r.dot( r0 ); - GEOSX_KRYLOV_BREAKDOWN_IF_ZERO( rho_old ) - GEOSX_KRYLOV_BREAKDOWN_IF_ZERO( omega ) + GEOS_KRYLOV_BREAKDOWN_IF_ZERO( rho_old ) + GEOS_KRYLOV_BREAKDOWN_IF_ZERO( omega ) // Compute beta real64 const beta = rho / rho_old * alpha / omega; @@ -107,7 +109,7 @@ void BicgstabSolver< VECTOR >::solve( Vector const & b, // Compute alpha real64 const vr0 = v.dot( r0 ); - GEOSX_KRYLOV_BREAKDOWN_IF_ZERO( vr0 ) + GEOS_KRYLOV_BREAKDOWN_IF_ZERO( vr0 ) alpha = rho / vr0; // compute x = x + alpha*y @@ -125,7 +127,7 @@ void BicgstabSolver< VECTOR >::solve( Vector const & b, // Update omega real64 const t2 = t.dot( t ); - GEOSX_KRYLOV_BREAKDOWN_IF_ZERO( t2 ) + GEOS_KRYLOV_BREAKDOWN_IF_ZERO( t2 ) omega = t.dot( s ) / t2; // Update x = x + omega*z diff --git a/src/coreComponents/linearAlgebra/solvers/BicgstabSolver.hpp b/src/coreComponents/linearAlgebra/solvers/BicgstabSolver.hpp index 0c940c7edfe..d1c9ed24d20 100644 --- a/src/coreComponents/linearAlgebra/solvers/BicgstabSolver.hpp +++ b/src/coreComponents/linearAlgebra/solvers/BicgstabSolver.hpp @@ -34,7 +34,7 @@ namespace geos * from Y. Saad (2003). */ template< typename VECTOR > -class BicgstabSolver : public KrylovSolver< VECTOR > +class BicgstabSolver final : public KrylovSolver< VECTOR > { public: @@ -71,9 +71,9 @@ class BicgstabSolver : public KrylovSolver< VECTOR > * @param [in] b system right hand side. * @param [inout] x system solution (input = initial guess, output = solution). */ - virtual void solve( Vector const & b, Vector & x ) const override final; + virtual void solve( Vector const & b, Vector & x ) const override; - virtual string methodName() const override final + virtual string methodName() const override { return "BiCGStab"; }; diff --git a/src/coreComponents/linearAlgebra/solvers/BlockPreconditioner.cpp b/src/coreComponents/linearAlgebra/solvers/BlockPreconditioner.cpp index 15840c09da6..4fa99cd1f4b 100644 --- a/src/coreComponents/linearAlgebra/solvers/BlockPreconditioner.cpp +++ b/src/coreComponents/linearAlgebra/solvers/BlockPreconditioner.cpp @@ -27,13 +27,9 @@ namespace geos { template< typename LAI > -BlockPreconditioner< LAI >::BlockPreconditioner( BlockShapeOption const shapeOption, - SchurComplementOption const schurOption, - BlockScalingOption const scalingOption ) +BlockPreconditioner< LAI >::BlockPreconditioner( LinearSolverParameters::Block params ) : Base(), - m_shapeOption( shapeOption ), - m_schurOption( schurOption ), - m_scalingOption( scalingOption ), + m_params( params ), m_matBlocks( 2, 2 ), m_solvers{}, m_scaling{ 1.0, 1.0 }, @@ -42,24 +38,26 @@ BlockPreconditioner< LAI >::BlockPreconditioner( BlockShapeOption const shapeOpt {} template< typename LAI > -BlockPreconditioner< LAI >::~BlockPreconditioner() = default; - -template< typename LAI > -void BlockPreconditioner< LAI >::reinitialize( Matrix const & mat, DofManager const & dofManager ) +void BlockPreconditioner< LAI >::reinitialize( Matrix const & mat ) { MPI_Comm const & comm = mat.comm(); if( m_blockDofs[1].empty() ) { - m_blockDofs[1] = dofManager.filterDofs( m_blockDofs[0] ); + GEOS_LAI_ASSERT( mat.dofManager() != nullptr ); + m_blockDofs[1] = mat.dofManager()->filterDofs( m_blockDofs[0] ); } for( localIndex i = 0; i < 2; ++i ) { - dofManager.makeRestrictor( m_blockDofs[i], comm, false, m_restrictors[i] ); - dofManager.makeRestrictor( m_blockDofs[i], comm, true, m_prolongators[i] ); - m_rhs( i ).create( m_restrictors[i].numLocalRows(), comm ); - m_sol( i ).create( m_restrictors[i].numLocalRows(), comm ); + if( m_prolongators[i] == nullptr ) + { + GEOS_LAI_ASSERT( mat.dofManager() != nullptr ); + mat.dofManager()->makeRestrictor( m_blockDofs[i], comm, true, m_prolongatorsOwned[i] ); + m_prolongators[i] = &m_prolongatorsOwned[i]; + } + m_rhs( i ).create( m_prolongators[i]->numLocalCols(), comm ); + m_sol( i ).create( m_prolongators[i]->numLocalCols(), comm ); } } @@ -75,16 +73,44 @@ void BlockPreconditioner< LAI >::setupBlock( localIndex const blockIndex, GEOS_LAI_ASSERT_GT( scaling, 0.0 ); m_blockDofs[blockIndex] = std::move( blockDofs ); - m_solvers[blockIndex] = std::move( solver ); + m_solversOwned[blockIndex] = std::move( solver ); + m_solvers[blockIndex] = m_solversOwned[blockIndex].get(); + m_scaling[blockIndex] = scaling; +} + +template< typename LAI > +void BlockPreconditioner< LAI >::setupBlock( localIndex const blockIndex, + std::vector< DofManager::SubComponent > blockDofs, + PreconditionerBase< LAI > * const solver, + real64 const scaling ) +{ + GEOS_LAI_ASSERT_GT( 2, blockIndex ); + GEOS_LAI_ASSERT( solver ); + GEOS_LAI_ASSERT( !blockDofs.empty() ); + GEOS_LAI_ASSERT_GT( scaling, 0.0 ); + + m_blockDofs[blockIndex] = std::move( blockDofs ); + m_solversOwned[blockIndex].reset(); + m_solvers[blockIndex] = solver; m_scaling[blockIndex] = scaling; } +template< typename LAI > +void BlockPreconditioner< LAI >::setProlongation( localIndex const blockIndex, + Matrix const & P ) +{ + GEOS_LAI_ASSERT_GT( 2, blockIndex ); + + m_prolongatorsOwned[blockIndex].reset(); + m_prolongators[blockIndex] = &P; +} + template< typename LAI > void BlockPreconditioner< LAI >::applyBlockScaling() { - if( m_scalingOption != BlockScalingOption::None ) + if( m_params.scaling != LinearSolverParameters::Block::Scaling::None ) { - if( m_scalingOption == BlockScalingOption::FrobeniusNorm ) + if( m_params.scaling == LinearSolverParameters::Block::Scaling::FrobeniusNorm ) { 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 ); @@ -104,14 +130,14 @@ void BlockPreconditioner< LAI >::applyBlockScaling() template< typename LAI > void BlockPreconditioner< LAI >::computeSchurComplement() { - switch( m_schurOption ) + switch( m_params.schurType ) { - case SchurComplementOption::None: + case LinearSolverParameters::Block::SchurType::None: { // nothing to do break; } - case SchurComplementOption::FirstBlockDiagonal: + case LinearSolverParameters::Block::SchurType::FirstBlockDiagonal: { m_matBlocks( 0, 0 ).extractDiagonal( m_rhs( 0 ) ); m_rhs( 0 ).reciprocal(); @@ -124,7 +150,7 @@ void BlockPreconditioner< LAI >::computeSchurComplement() m_matBlocks( 0, 1 ).leftScale( m_rhs( 0 ) ); break; } - case SchurComplementOption::RowsumDiagonalProbing: + case LinearSolverParameters::Block::SchurType::RowsumDiagonalProbing: { m_sol( 1 ).set( -1.0 ); m_matBlocks( 0, 1 ).apply( m_sol( 1 ), m_rhs( 0 ) ); @@ -133,7 +159,7 @@ void BlockPreconditioner< LAI >::computeSchurComplement() m_matBlocks( 1, 1 ).addDiagonal( m_rhs( 1 ), 1.0 ); break; } - case SchurComplementOption::FirstBlockUserDefined: + case LinearSolverParameters::Block::SchurType::FirstBlockUserDefined: { Matrix const & prec00 = m_solvers[0]->preconditionerMatrix(); Matrix mat11; @@ -151,9 +177,6 @@ void BlockPreconditioner< LAI >::computeSchurComplement() template< typename LAI > void BlockPreconditioner< LAI >::setup( Matrix const & mat ) { - // Check that DofManager is available - GEOS_LAI_ASSERT_MSG( mat.dofManager() != nullptr, "BlockPreconditioner requires a DofManager" ); - // Check that user has set block solvers GEOS_LAI_ASSERT( m_solvers[0] != nullptr ); GEOS_LAI_ASSERT( m_solvers[1] != nullptr ); @@ -170,18 +193,25 @@ void BlockPreconditioner< LAI >::setup( Matrix const & mat ) // If the matrix size/structure has changed, need to resize internal LA objects and recompute restrictors. if( newSize ) { - reinitialize( mat, *mat.dofManager() ); + reinitialize( mat ); } // Extract diagonal blocks - mat.multiplyPtAP( m_prolongators[0], m_matBlocks( 0, 0 ) ); - mat.multiplyPtAP( m_prolongators[1], m_matBlocks( 1, 1 ) ); + mat.multiplyPtAP( *m_prolongators[0], m_matBlocks( 0, 0 ) ); + mat.multiplyPtAP( *m_prolongators[1], m_matBlocks( 1, 1 ) ); + + // HACK: a coupled DofManager is technically not compatible with diagonal blocks. + // We can create "reduced" managers using DofManager::setupFrom(), but it can be a waste of time, + // Instead, we expect block solvers to not rely on global information and instead query by field. + m_matBlocks( 0, 0 ).setDofManager( mat.dofManager() ); + m_matBlocks( 1, 1 ).setDofManager( mat.dofManager() ); // Extract off-diagonal blocks only if used - if( m_schurOption != SchurComplementOption::None && m_shapeOption != BlockShapeOption::Diagonal ) + if( m_params.schurType != LinearSolverParameters::Block::SchurType::None + || m_params.shape != LinearSolverParameters::Block::Shape::Diagonal ) { - mat.multiplyRAP( m_restrictors[0], m_prolongators[1], m_matBlocks( 0, 1 ) ); - mat.multiplyRAP( m_restrictors[1], m_prolongators[0], m_matBlocks( 1, 0 ) ); + mat.multiplyPtAP( *m_prolongators[0], *m_prolongators[1], m_matBlocks( 0, 1 ) ); + mat.multiplyPtAP( *m_prolongators[1], *m_prolongators[0], m_matBlocks( 1, 0 ) ); } applyBlockScaling(); @@ -194,36 +224,41 @@ template< typename LAI > void BlockPreconditioner< LAI >::apply( Vector const & src, Vector & dst ) const { - m_restrictors[0].apply( src, m_rhs( 0 ) ); - m_restrictors[1].apply( src, m_rhs( 1 ) ); + using Shape = LinearSolverParameters::Block::Shape; + + m_prolongators[0]->applyTranspose( src, m_rhs( 0 ) ); + m_prolongators[1]->applyTranspose( src, m_rhs( 1 ) ); for( localIndex i = 0; i < 2; ++i ) { m_rhs( i ).scale( m_scaling[i] ); } - // Perform a predictor step by solving (0,0) block and subtracting from 1-block rhs - if( m_shapeOption == BlockShapeOption::LowerUpperTriangular ) + if( m_params.shape == Shape::LowerUpperTriangular || m_params.shape == Shape::LowerTriangular ) { + // Solve the (0,0)-block and update (1,1)-block rhs m_solvers[0]->apply( m_rhs( 0 ), m_sol( 0 ) ); m_matBlocks( 1, 0 ).residual( m_sol( 0 ), m_rhs( 1 ), m_rhs( 1 ) ); } - // Solve the (1,1) block modified via Schur complement + // Solve the (1,1)-block modified via Schur complement m_solvers[1]->apply( m_rhs( 1 ), m_sol( 1 ) ); - // Update the 0-block rhs - if( m_shapeOption != BlockShapeOption::Diagonal ) + if( m_params.shape != Shape::LowerTriangular ) { - m_matBlocks( 0, 1 ).residual( m_sol( 1 ), m_rhs( 0 ), m_rhs( 0 ) ); - } + if( m_params.shape != Shape::Diagonal ) + { + // Update the 0-block rhs + m_matBlocks( 0, 1 ).residual( m_sol( 1 ), m_rhs( 0 ), m_rhs( 0 ) ); + } - // Solve the (0,0) block with the current rhs - m_solvers[0]->apply( m_rhs( 0 ), m_sol( 0 ) ); + // Solve the (0,0) block with the current rhs + m_solvers[0]->apply( m_rhs( 0 ), m_sol( 0 ) ); + } // Combine block solutions into global solution vector - m_prolongators[0].apply( m_sol( 0 ), dst ); - m_prolongators[1].gemv( 1.0, m_sol( 1 ), 1.0, dst ); + m_prolongators[0]->apply( m_sol( 0 ), dst ); + m_prolongators[1]->gemv( 1.0, m_sol( 1 ), 1.0, dst ); } template< typename LAI > @@ -232,8 +267,8 @@ void BlockPreconditioner< LAI >::clear() Base::clear(); for( localIndex i = 0; i < 2; ++i ) { - m_restrictors[i].reset(); - m_prolongators[i].reset(); + m_prolongatorsOwned[i].reset(); + m_prolongators[i] = nullptr; m_solvers[i]->clear(); m_rhs( i ).reset(); m_sol( i ).reset(); diff --git a/src/coreComponents/linearAlgebra/solvers/BlockPreconditioner.hpp b/src/coreComponents/linearAlgebra/solvers/BlockPreconditioner.hpp index f6ea90dd635..f0f90545188 100644 --- a/src/coreComponents/linearAlgebra/solvers/BlockPreconditioner.hpp +++ b/src/coreComponents/linearAlgebra/solvers/BlockPreconditioner.hpp @@ -23,43 +23,11 @@ #include "linearAlgebra/common/PreconditionerBase.hpp" #include "linearAlgebra/utilities/BlockOperator.hpp" #include "linearAlgebra/utilities/BlockVector.hpp" +#include "linearAlgebra/utilities/LinearSolverParameters.hpp" namespace geos { -/** - * @brief Type of Schur complement approximation used - * - * @todo Need more descriptive names for options - */ -enum class SchurComplementOption -{ - None, //!< No Schur complement - just block-GS/block-Jacobi preconditioner - FirstBlockDiagonal, //!< Approximate first block with its diagonal - RowsumDiagonalProbing, //!< Rowsum-preserving diagonal approximation constructed with probing - FirstBlockUserDefined //!< User defined preconditioner for the first block -}; - -/** - * @brief Type of block row scaling to apply - */ -enum class BlockScalingOption -{ - None, //!< No scaling - FrobeniusNorm, //!< Equilibrate Frobenius norm of the diagonal blocks - UserProvided //!< User-provided scaling -}; - -/** - * @brief Shape of the block preconditioner - */ -enum class BlockShapeOption -{ - Diagonal, //!< (D)^{-1} - UpperTriangular, //!< (DU)^{-1} - LowerUpperTriangular //!< (LDU)^{-1} -}; - /* * Since formulas in Doxygen are broken with 1.8.13 and certain versions of ghostscript, * keeping this documentation in a separate comment block for now. Should be moved into @@ -112,18 +80,25 @@ class BlockPreconditioner : public PreconditionerBase< LAI > /** * @brief Constructor. - * @param shapeOption preconditioner block shape - * @param schurOption type of Schur complement approximation to use - * @param scalingOption type of scaling to apply to blocks + * @param params block preconditioner parameters */ - explicit BlockPreconditioner( BlockShapeOption const shapeOption, - SchurComplementOption const schurOption, - BlockScalingOption const scalingOption ); + explicit BlockPreconditioner( LinearSolverParameters::Block params ); /** - * @brief Destructor. + * @brief Setup data for one of the two blocks. + * @param blockIndex index of the block to set up + * @param blockDofs choice of DoF components (from a monolithic system) + * @param solver instance of the inner preconditioner for the block + * @param scaling user-provided row scaling coefficient for this block + * + * @note While not strictly required, it is generally expected that @p blockDofs + * of the two blocks are non-overlapping and their union includes all + * DoF components in the monolithic system. */ - virtual ~BlockPreconditioner() override; + void setupBlock( localIndex const blockIndex, + std::vector< DofManager::SubComponent > blockDofs, + std::unique_ptr< PreconditionerBase< LAI > > solver, + real64 const scaling = 1.0 ); /** * @brief Setup data for one of the two blocks. @@ -132,22 +107,30 @@ class BlockPreconditioner : public PreconditionerBase< LAI > * @param solver instance of the inner preconditioner for the block * @param scaling user-provided row scaling coefficient for this block * + * @note This version does not take ownership of @p solver that must have its lifetime guaranteed elsewhere. + * * @note While not strictly required, it is generally expected that @p blockDofs * of the two blocks are non-overlapping and their union includes all * DoF components in the monolithic system. */ void setupBlock( localIndex const blockIndex, std::vector< DofManager::SubComponent > blockDofs, - std::unique_ptr< PreconditionerBase< LAI > > solver, + PreconditionerBase< LAI > * const solver, real64 const scaling = 1.0 ); + /** + * @brief Set user-provided prolongation operator for a block + * @param blockIndex index of the block to set up + * @param P the operator + */ + void setProlongation( localIndex const blockIndex, + Matrix const & P ); + /** * @name PreconditionerBase interface methods */ ///@{ - using PreconditionerBase< LAI >::setup; - /** * @brief Compute the preconditioner from a matrix * @param mat the matrix to precondition @@ -167,14 +150,21 @@ class BlockPreconditioner : public PreconditionerBase< LAI > ///@} + /** + * @brief @return the block operator consisting of extracted matrix blocks + */ + BlockOperator< Vector, Matrix > const & blocks() const + { + return m_matBlocks; + } + private: /** * @brief Initialize/resize internal data structures for a new linear system. * @param mat the new system matrix - * @param dofManager the new dof manager */ - void reinitialize( Matrix const & mat, DofManager const & dofManager ); + void reinitialize( Matrix const & mat ); /** * @brief Apply block scaling to system blocks (which must be already extracted). @@ -186,32 +176,29 @@ class BlockPreconditioner : public PreconditionerBase< LAI > */ void computeSchurComplement(); - /// Shape of the block preconditioner - BlockShapeOption m_shapeOption; - - /// Type of Schur complement to construct - SchurComplementOption m_schurOption; - - /// Whether to scale blocks to equilibrate norms - BlockScalingOption m_scalingOption; + /// Block preconditioner parameters + LinearSolverParameters::Block m_params; /// Description of dof components making up each of the two main blocks - std::array< std::vector< DofManager::SubComponent >, 2 > m_blockDofs; + std::array< std::vector< DofManager::SubComponent >, 2 > m_blockDofs{}; - /// Restriction operators for each sub-block - std::array< Matrix, 2 > m_restrictors; + /// Pointers to prolongation operators for each sub-block + std::array< Matrix const *, 2 > m_prolongators{}; /// Prolongation operators for each sub-block - std::array< Matrix, 2 > m_prolongators; + std::array< Matrix, 2 > m_prolongatorsOwned{}; /// Matrix blocks BlockOperator< Vector, Matrix > m_matBlocks; - /// Individual block preconditioners - std::array< std::unique_ptr< PreconditionerBase< LAI > >, 2 > m_solvers; + /// Individual block preconditioner pointers + std::array< PreconditionerBase< LAI > *, 2 > m_solvers{}; + + /// Individual block preconditioner operators (when owned) + std::array< std::unique_ptr< PreconditionerBase< LAI > >, 2 > m_solversOwned{}; /// Scaling of each block - std::array< real64, 2 > m_scaling; + std::array< real64, 2 > m_scaling{}; /// Internal vector of block residuals mutable BlockVector< Vector > m_rhs; diff --git a/src/coreComponents/linearAlgebra/solvers/CgSolver.cpp b/src/coreComponents/linearAlgebra/solvers/CgSolver.cpp index bceaa44c1da..9ffc00d18f3 100644 --- a/src/coreComponents/linearAlgebra/solvers/CgSolver.cpp +++ b/src/coreComponents/linearAlgebra/solvers/CgSolver.cpp @@ -20,6 +20,7 @@ #include "CgSolver.hpp" +#include "common/TimingMacros.hpp" #include "common/Stopwatch.hpp" #include "linearAlgebra/interfaces/InterfaceTypes.hpp" #include "common/LinearOperator.hpp" @@ -57,6 +58,7 @@ CgSolver< VECTOR >::CgSolver( LinearSolverParameters params, template< typename VECTOR > void CgSolver< VECTOR >::solve( Vector const & b, Vector & x ) const { + GEOS_MARK_FUNCTION; Stopwatch watch; // Define residual vector @@ -113,7 +115,7 @@ void CgSolver< VECTOR >::solve( Vector const & b, Vector & x ) const // compute alpha real64 const pAp = p.dot( Ap ); - GEOSX_KRYLOV_BREAKDOWN_IF_ZERO( pAp ) + GEOS_KRYLOV_BREAKDOWN_IF_ZERO( pAp ) real64 const alpha = tau / pAp; // Update x = x + alpha*p diff --git a/src/coreComponents/linearAlgebra/solvers/CgSolver.hpp b/src/coreComponents/linearAlgebra/solvers/CgSolver.hpp index f5d6b821bb0..540a29351f3 100644 --- a/src/coreComponents/linearAlgebra/solvers/CgSolver.hpp +++ b/src/coreComponents/linearAlgebra/solvers/CgSolver.hpp @@ -34,7 +34,7 @@ namespace geos * from Y. Saad (2003). */ template< typename VECTOR > -class CgSolver : public KrylovSolver< VECTOR > +class CgSolver final : public KrylovSolver< VECTOR > { public: @@ -71,9 +71,9 @@ class CgSolver : public KrylovSolver< VECTOR > * @param [in] b system right hand side. * @param [inout] x system solution (input = initial guess, output = solution). */ - virtual void solve( Vector const & b, Vector & x ) const override final; + virtual void solve( Vector const & b, Vector & x ) const override; - virtual string methodName() const override final + virtual string methodName() const override { return "CG"; }; @@ -96,6 +96,6 @@ class CgSolver : public KrylovSolver< VECTOR > }; -} // namespace GEOSX +} // namespace geos #endif /*GEOS_LINEARALGEBRA_SOLVERS_CGSOLVER_HPP_*/ diff --git a/src/coreComponents/linearAlgebra/solvers/GmresSolver.cpp b/src/coreComponents/linearAlgebra/solvers/GmresSolver.cpp index 4274c50f294..4d3620a93c8 100644 --- a/src/coreComponents/linearAlgebra/solvers/GmresSolver.cpp +++ b/src/coreComponents/linearAlgebra/solvers/GmresSolver.cpp @@ -18,6 +18,7 @@ #include "GmresSolver.hpp" +#include "common/TimingMacros.hpp" #include "common/Stopwatch.hpp" #include "linearAlgebra/interfaces/InterfaceTypes.hpp" #include "linearAlgebra/solvers/KrylovUtils.hpp" @@ -30,11 +31,10 @@ template< typename VECTOR > GmresSolver< VECTOR >::GmresSolver( LinearSolverParameters params, LinearOperator< Vector > const & A, LinearOperator< Vector > const & M ) - : KrylovSolver< VECTOR >( std::move( params ), A, M ), - m_kspace( m_params.krylov.maxRestart + 1 ), - m_kspaceInitialized( false ) + : KrylovSolver< VECTOR >( std::move( params ), A, M ) { GEOS_ERROR_IF_LE_MSG( m_params.krylov.maxRestart, 0, "GMRES: max number of iterations until restart must be positive." ); + m_kspace.reserve( m_params.krylov.maxRestart + 1 ); } namespace @@ -89,16 +89,21 @@ template< typename VECTOR > void GmresSolver< VECTOR >::solve( Vector const & b, Vector & x ) const { - // We create Krylov subspace vectors once using the size and partitioning of b. - // On repeated calls to solve() input vectors must have the same size and partitioning. - if( !m_kspaceInitialized ) + GEOS_MARK_FUNCTION; + + // Function to grow K-space on demand + auto const addKspaceVector = [&]( integer const n ) { - for( VectorTemp & kv : m_kspace ) + for( integer i = m_kspace.size(); i <= n; ++i ) { - kv = createTempVector( b ); + m_kspace.emplace_back( createTempVector( b ) ); } - m_kspaceInitialized = true; - } + }; + + // TEMP: evaluate the speed benefit of pre-allocation vs memory waste. + // Rationale: LA packages (e.g. hypre) allocate krylov space this during + // setup, so we don't want to include it in the timing of the solve. + addKspaceVector( m_params.krylov.maxRestart ); Stopwatch watch; @@ -126,6 +131,10 @@ void GmresSolver< VECTOR >::solve( Vector const & b, m_result.status = LinearSolverResult::Status::NotConverged; m_residualNorms.clear(); + // On subsequent calls to solve(), b must have the same size/distribution + addKspaceVector( 0 ); + GEOS_LAI_ASSERT_EQ( b.localSize(), m_kspace[0].localSize() ); + integer & k = m_result.numIterations; while( k <= m_params.krylov.maxIterations && m_result.status == LinearSolverResult::Status::NotConverged ) { @@ -165,7 +174,9 @@ void GmresSolver< VECTOR >::solve( Vector const & b, } H( j+1, j ) = w.norm2(); - GEOSX_KRYLOV_BREAKDOWN_IF_ZERO( H( j+1, j ) ) + GEOS_KRYLOV_BREAKDOWN_IF_ZERO( H( j+1, j ) ) + + addKspaceVector( j+1 ); m_kspace[j+1].axpby( 1.0 / H( j+1, j ), w, 0.0 ); // Apply all previous rotations to the new column diff --git a/src/coreComponents/linearAlgebra/solvers/GmresSolver.hpp b/src/coreComponents/linearAlgebra/solvers/GmresSolver.hpp index 15edabbe0a4..2dea7e9f32a 100644 --- a/src/coreComponents/linearAlgebra/solvers/GmresSolver.hpp +++ b/src/coreComponents/linearAlgebra/solvers/GmresSolver.hpp @@ -34,7 +34,7 @@ namespace geos * from Y. Saad (2003). */ template< typename VECTOR > -class GmresSolver : public KrylovSolver< VECTOR > +class GmresSolver final : public KrylovSolver< VECTOR > { public: @@ -71,9 +71,9 @@ class GmresSolver : public KrylovSolver< VECTOR > * @param [in] b system right hand side. * @param [inout] x system solution (input = initial guess, output = solution). */ - virtual void solve( Vector const & b, Vector & x ) const override final; + virtual void solve( Vector const & b, Vector & x ) const override; - virtual string methodName() const override final + virtual string methodName() const override { return "GMRES"; }; @@ -95,10 +95,7 @@ class GmresSolver : public KrylovSolver< VECTOR > using Base::logResult; /// Storage for Krylov subspace vectors - array1d< VectorTemp > m_kspace; - - /// Flag indicating whether kspace vectors have been created - bool mutable m_kspaceInitialized; + mutable array1d< VectorTemp > m_kspace; }; } // namespace geos diff --git a/src/coreComponents/linearAlgebra/solvers/KrylovSolver.cpp b/src/coreComponents/linearAlgebra/solvers/KrylovSolver.cpp index d3c461ba08b..2148df1f641 100644 --- a/src/coreComponents/linearAlgebra/solvers/KrylovSolver.cpp +++ b/src/coreComponents/linearAlgebra/solvers/KrylovSolver.cpp @@ -20,6 +20,7 @@ #include "linearAlgebra/solvers/BicgstabSolver.hpp" #include "linearAlgebra/solvers/CgSolver.hpp" #include "linearAlgebra/solvers/GmresSolver.hpp" +#include "linearAlgebra/solvers/RichardsonSolver.hpp" #include "linearAlgebra/interfaces/InterfaceTypes.hpp" namespace geos @@ -67,6 +68,12 @@ KrylovSolver< VECTOR >::create( LinearSolverParameters const & parameters, matrix, precond ); } + case LinearSolverParameters::SolverType::richardson: + { + return std::make_unique< RichardsonSolver< Vector > >( parameters, + matrix, + precond ); + } default: { GEOS_ERROR( "Unsupported linear solver type: " << parameters.solverType ); @@ -81,8 +88,18 @@ void KrylovSolver< VECTOR >::logProgress() const GEOS_ASSERT( !m_residualNorms.empty() ); if( m_params.logLevel >= 2 ) { - real64 const relNorm = m_residualNorms[0] > 0.0 ? m_residualNorms.back() / m_residualNorms[0] : 0.0; - GEOS_LOG_RANK_0( GEOS_FMT( "[{}] iteration {}: residual = {:e}", methodName(), m_result.numIterations, relNorm ) ); + constexpr char const * headFormat = "{:>4} {:>12} {:>9} {:>12}"; + constexpr char const * lineFormat = "{:>4} {:>12.6e} {:>9.6f} {:>12.6e}"; + integer const iter = m_result.numIterations; + if( iter == 0 ) + { + GEOS_LOG_RANK_0( GEOS_FMT( "[{}] start iteration", methodName() ) ); + GEOS_LOG_RANK_0( GEOS_FMT( headFormat, "iter", "resid.norm", "conv.rate", "rel.res.norm" ) ); + } + real64 const norm = m_residualNorms[iter]; + real64 const relNorm = m_residualNorms[0] > 0.0 ? norm / m_residualNorms[0] : 0.0; + real64 const convRate = iter > 0 ? norm / m_residualNorms[iter - 1] : 1.0; + GEOS_LOG_RANK_0( GEOS_FMT( lineFormat, iter, norm, convRate, relNorm ) ); } } diff --git a/src/coreComponents/linearAlgebra/solvers/KrylovSolver.hpp b/src/coreComponents/linearAlgebra/solvers/KrylovSolver.hpp index 33504c4b105..bab8c9932a5 100644 --- a/src/coreComponents/linearAlgebra/solvers/KrylovSolver.hpp +++ b/src/coreComponents/linearAlgebra/solvers/KrylovSolver.hpp @@ -12,6 +12,10 @@ * ------------------------------------------------------------------------------------------------------------ */ +/** + * @file KrylovSolver.hpp + */ + #ifndef GEOS_LINEARALGEBRA_SOLVERS_KRYLOVSOLVER_HPP_ #define GEOS_LINEARALGEBRA_SOLVERS_KRYLOVSOLVER_HPP_ @@ -46,9 +50,10 @@ class KrylovSolver : public LinearOperator< VECTOR > * @param precond preconditioning operator (must be set up by the user prior to calling solve()/apply()) * @return an owning pointer to the newly instantiated solver */ - static std::unique_ptr< KrylovSolver< VECTOR > > create( LinearSolverParameters const & parameters, - LinearOperator< VECTOR > const & matrix, - LinearOperator< VECTOR > const & precond ); + static std::unique_ptr< KrylovSolver< VECTOR > > + create( LinearSolverParameters const & parameters, + LinearOperator< VECTOR > const & matrix, + LinearOperator< VECTOR > const & precond ); /** * @brief Constructor. @@ -60,11 +65,6 @@ class KrylovSolver : public LinearOperator< VECTOR > LinearOperator< Vector > const & matrix, LinearOperator< Vector > const & precond ); - /** - * @brief Virtual destructor - */ - virtual ~KrylovSolver() override = default; - /** * @brief Solve preconditioned system * @param [in] b system right hand side. @@ -72,10 +72,8 @@ class KrylovSolver : public LinearOperator< VECTOR > */ virtual void solve( Vector const & b, Vector & x ) const = 0; - /** * @brief Apply operator to a vector. - * * @param src Input vector (src). * @param dst Output vector (dst). */ diff --git a/src/coreComponents/linearAlgebra/solvers/KrylovUtils.hpp b/src/coreComponents/linearAlgebra/solvers/KrylovUtils.hpp index 01bf8aef804..8a7a34427a1 100644 --- a/src/coreComponents/linearAlgebra/solvers/KrylovUtils.hpp +++ b/src/coreComponents/linearAlgebra/solvers/KrylovUtils.hpp @@ -24,7 +24,7 @@ * @brief Exit solver iteration and report a breakdown if value too close to zero. * @param VAR the variable or expression */ -#define GEOSX_KRYLOV_BREAKDOWN_IF_ZERO( VAR ) \ +#define GEOS_KRYLOV_BREAKDOWN_IF_ZERO( VAR ) \ if( isZero( VAR, 0.0 ) ) \ { \ if( m_params.logLevel >= 1 ) \ diff --git a/src/coreComponents/linearAlgebra/solvers/PreconditionerIdentity.hpp b/src/coreComponents/linearAlgebra/solvers/PreconditionerIdentity.hpp index 1d3a9e133c4..a48b42d81ce 100644 --- a/src/coreComponents/linearAlgebra/solvers/PreconditionerIdentity.hpp +++ b/src/coreComponents/linearAlgebra/solvers/PreconditionerIdentity.hpp @@ -12,10 +12,13 @@ * ------------------------------------------------------------------------------------------------------------ */ -#ifndef GEOS_LINEARALGEBRA_SOLVERS_PRECONDITIONERIDENTITY_HPP_ -#define GEOS_LINEARALGEBRA_SOLVERS_PRECONDITIONERIDENTITY_HPP_ +/** + * @file PreconditionerIdentity.hpp + */ + +#ifndef GEOS_LINEARALGEBRA_SOLVERS_PRECONDITIONERNULL_HPP_ +#define GEOS_LINEARALGEBRA_SOLVERS_PRECONDITIONERNULL_HPP_ -#include "linearAlgebra/common/LinearOperator.hpp" #include "linearAlgebra/common/PreconditionerBase.hpp" namespace geos @@ -39,8 +42,6 @@ class PreconditionerIdentity : public PreconditionerBase< LAI > /// Alias for matrix type using Matrix = typename Base::Matrix; - virtual ~PreconditionerIdentity() = default; - /** * @brief Apply operator to a vector. * @@ -56,6 +57,6 @@ class PreconditionerIdentity : public PreconditionerBase< LAI > } }; -} +} // namespace geos #endif //GEOS_LINEARALGEBRA_SOLVERS_PRECONDITIONERIDENTITY_HPP_ diff --git a/src/coreComponents/linearAlgebra/solvers/PreconditionerNull.hpp b/src/coreComponents/linearAlgebra/solvers/PreconditionerNull.hpp new file mode 100644 index 00000000000..21389a738ad --- /dev/null +++ b/src/coreComponents/linearAlgebra/solvers/PreconditionerNull.hpp @@ -0,0 +1,62 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file PreconditionerNull.hpp + */ + +#ifndef GEOSX_LINEARALGEBRA_SOLVERS_PRECONDITIONERNULL_HPP_ +#define GEOSX_LINEARALGEBRA_SOLVERS_PRECONDITIONERNULL_HPP_ + +#include "linearAlgebra/common/PreconditionerBase.hpp" + +namespace geos +{ + +/** + * @brief Common interface for identity preconditioning operator + * @tparam LAI linear algebra interface providing vectors, matrices and solvers + */ +template< typename LAI > +class PreconditionerNull : public PreconditionerBase< LAI > +{ +public: + + /// Alias for base type + using Base = PreconditionerBase< LAI >; + + /// Alias for vector type + using Vector = typename Base::Vector; + + /// Alias for matrix type + using Matrix = typename Base::Matrix; + + /** + * @brief Apply operator to a vector. + * + * @param src Input vector (src). + * @param dst Output vector (dst). + */ + virtual void apply( Vector const & src, + Vector & dst ) const override + { + GEOS_LAI_ASSERT_EQ( this->numGlobalRows(), dst.globalSize() ); + GEOS_LAI_ASSERT_EQ( this->numGlobalCols(), src.globalSize() ); + dst.zero(); + } +}; + +} // namespace geos + +#endif //GEOSX_LINEARALGEBRA_SOLVERS_PRECONDITIONERNULL_HPP_ diff --git a/src/coreComponents/linearAlgebra/solvers/RichardsonSolver.cpp b/src/coreComponents/linearAlgebra/solvers/RichardsonSolver.cpp new file mode 100644 index 00000000000..aea01f3b4f3 --- /dev/null +++ b/src/coreComponents/linearAlgebra/solvers/RichardsonSolver.cpp @@ -0,0 +1,100 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file RichardsonSolver.cpp + */ + +#include "RichardsonSolver.hpp" + +#include "common/TimingMacros.hpp" +#include "common/Stopwatch.hpp" +#include "linearAlgebra/interfaces/InterfaceTypes.hpp" + +namespace geos +{ + +template< typename VECTOR > +RichardsonSolver< VECTOR >::RichardsonSolver( LinearSolverParameters params, + LinearOperator< Vector > const & A, + LinearOperator< Vector > const & M ) + : KrylovSolver< VECTOR >( std::move( params ), A, M ), + m_omega( m_params.relaxation.weight ) +{} + +template< typename VECTOR > +void RichardsonSolver< VECTOR >::solve( Vector const & b, + Vector & x ) const +{ + GEOS_MARK_FUNCTION; + Stopwatch watch; + + VectorTemp r = createTempVector( b ); + VectorTemp z = createTempVector( b ); + + // Compute initial rk + m_operator.residual( x, b, r ); + + // Compute the target absolute tolerance + real64 const rnorm0 = r.norm2(); + real64 const absTol = rnorm0 * m_params.krylov.relTolerance; + + // Initialize iteration state + m_result.status = LinearSolverResult::Status::NotConverged; + m_residualNorms.clear(); + + integer & k = m_result.numIterations; + for( k = 0; k <= m_params.krylov.maxIterations; ++k ) + { + real64 const rnorm = r.norm2(); + m_residualNorms.emplace_back( rnorm ); + logProgress(); + + // Convergence check on ||rk||/||r0|| + if( rnorm <= absTol ) + { + m_result.status = LinearSolverResult::Status::Success; + break; + } + + // Update: z = Mr, x = x + w*z, r = b - Ax + m_precond.apply( r, z ); + x.axpy( m_omega, z ); + m_operator.residual( x, b, r ); + } + + m_result.residualReduction = rnorm0 > 0.0 ? m_residualNorms.back() / rnorm0 : 0.0; + m_result.solveTime = watch.elapsedTime(); + logResult(); +} + +// ----------------------- +// Explicit Instantiations +// ----------------------- +#ifdef GEOSX_USE_TRILINOS +template class RichardsonSolver< TrilinosInterface::ParallelVector >; +template class RichardsonSolver< BlockVectorView< TrilinosInterface::ParallelVector > >; +#endif + +#ifdef GEOSX_USE_HYPRE +template class RichardsonSolver< HypreInterface::ParallelVector >; +template class RichardsonSolver< BlockVectorView< HypreInterface::ParallelVector > >; +#endif + +#ifdef GEOSX_USE_PETSC +template class RichardsonSolver< PetscInterface::ParallelVector >; +template class RichardsonSolver< BlockVectorView< PetscInterface::ParallelVector > >; +#endif + +} // geosx diff --git a/src/coreComponents/linearAlgebra/solvers/RichardsonSolver.hpp b/src/coreComponents/linearAlgebra/solvers/RichardsonSolver.hpp new file mode 100644 index 00000000000..d483ca300dd --- /dev/null +++ b/src/coreComponents/linearAlgebra/solvers/RichardsonSolver.hpp @@ -0,0 +1,102 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file RichardsonSolver.hpp + */ + +#ifndef GEOSX_LINEARALGEBRA_SOLVERS_RICHARDSONSOLVER_HPP_ +#define GEOSX_LINEARALGEBRA_SOLVERS_RICHARDSONSOLVER_HPP_ + +#include "linearAlgebra/solvers/KrylovSolver.hpp" + +namespace geos +{ + +/** + * @brief Implements right-preconditioned modified Richardson iteration. + * @tparam VECTOR type of vectors this solver operates on + * @note Richardson is not a Krylov subspace method, but + * for convenience inherits from KrylovSolver; that + * class should really be renamed to IterativeSolver. + */ +template< typename VECTOR > +class RichardsonSolver final : public KrylovSolver< VECTOR > +{ +public: + + /// Alias for the base type + using Base = KrylovSolver< VECTOR >; + + /// Alias for the vector type + using Vector = typename Base::Vector; + + /** + * @name Constructor/Destructor Methods + */ + ///@{ + + /** + * @brief Solver object constructor. + * @param[in] params parameters for the solver + * @param[in] matrix reference to the system matrix + * @param[in] precond reference to the preconditioning operator + */ + RichardsonSolver( LinearSolverParameters params, + LinearOperator< Vector > const & matrix, + LinearOperator< Vector > const & precond ); + + ///@} + + /** + * @name KrylovSolver interface + */ + ///@{ + + /** + * @brief Solve preconditioned system + * @param [in] b system right hand side. + * @param [inout] x system solution (input = initial guess, output = solution). + */ + virtual void solve( Vector const & b, Vector & x ) const override; + + virtual string methodName() const override + { + return "Richardson"; + }; + + ///@} + +protected: + + /// Alias for vector type that can be used for temporaries + using VectorTemp = typename KrylovSolver< VECTOR >::VectorTemp; + + using Base::m_params; + using Base::m_operator; + using Base::m_precond; + using Base::m_residualNorms; + using Base::m_result; + using Base::createTempVector; + using Base::logProgress; + using Base::logResult; + +private: + + real64 m_omega; +}; + +} // geosx + +#endif //GEOSX_LINEARALGEBRA_SOLVERS_RICHARDSONSOLVER_HPP_ diff --git a/src/coreComponents/linearAlgebra/unitTests/testExternalSolvers.cpp b/src/coreComponents/linearAlgebra/unitTests/testExternalSolvers.cpp index 41e78945189..314d6d0a4d5 100644 --- a/src/coreComponents/linearAlgebra/unitTests/testExternalSolvers.cpp +++ b/src/coreComponents/linearAlgebra/unitTests/testExternalSolvers.cpp @@ -47,7 +47,7 @@ LinearSolverParameters params_GMRES_ILU() parameters.krylov.relTolerance = 1e-8; parameters.krylov.maxIterations = 300; parameters.solverType = LinearSolverParameters::SolverType::gmres; - parameters.preconditionerType = LinearSolverParameters::PreconditionerType::iluk; + parameters.preconditionerType = LinearSolverParameters::PreconditionerType::ilu; parameters.ifact.fill = 1; return parameters; } diff --git a/src/coreComponents/linearAlgebra/unitTests/testKrylovSolvers.cpp b/src/coreComponents/linearAlgebra/unitTests/testKrylovSolvers.cpp index da8d734d661..9efd8a3cc91 100644 --- a/src/coreComponents/linearAlgebra/unitTests/testKrylovSolvers.cpp +++ b/src/coreComponents/linearAlgebra/unitTests/testKrylovSolvers.cpp @@ -54,6 +54,16 @@ LinearSolverParameters params_GMRES() return parameters; } +LinearSolverParameters params_Richardson() +{ + LinearSolverParameters parameters; + parameters.krylov.relTolerance = 1e-4; + parameters.krylov.maxIterations = 500; + parameters.solverType = geos::LinearSolverParameters::SolverType::richardson; + parameters.relaxation.weight = 0.2; + return parameters; +} + template< typename OPERATOR, typename PRECOND, typename VECTOR > class KrylovSolverTestBase : public ::testing::Test { @@ -153,10 +163,16 @@ TYPED_TEST_P( KrylovSolverTest, GMRES ) this->test( params_GMRES() ); } +TYPED_TEST_P( KrylovSolverTest, Richardson ) +{ + this->test( params_Richardson() ); +} + REGISTER_TYPED_TEST_SUITE_P( KrylovSolverTest, CG, BiCGSTAB, - GMRES ); + GMRES, + Richardson ); #ifdef GEOSX_USE_TRILINOS INSTANTIATE_TYPED_TEST_SUITE_P( Trilinos, KrylovSolverTest, TrilinosInterface, ); @@ -243,10 +259,17 @@ TYPED_TEST_P( KrylovSolverBlockTest, GMRES ) this->test( params_GMRES() ); } +TYPED_TEST_P( KrylovSolverBlockTest, Richardson ) +{ + this->test( params_Richardson() ); +} + + REGISTER_TYPED_TEST_SUITE_P( KrylovSolverBlockTest, CG, BiCGSTAB, - GMRES ); + GMRES, + Richardson ); #ifdef GEOSX_USE_TRILINOS INSTANTIATE_TYPED_TEST_SUITE_P( Trilinos, KrylovSolverBlockTest, TrilinosInterface, ); diff --git a/src/coreComponents/linearAlgebra/unitTests/testLinearSolverParametersEnums.cpp b/src/coreComponents/linearAlgebra/unitTests/testLinearSolverParametersEnums.cpp index c136a58229c..97140425663 100644 --- a/src/coreComponents/linearAlgebra/unitTests/testLinearSolverParametersEnums.cpp +++ b/src/coreComponents/linearAlgebra/unitTests/testLinearSolverParametersEnums.cpp @@ -33,6 +33,7 @@ TEST( LinearSolverParametersEnums, SolverType ) ASSERT_EQ( "gmres", toString( EnumType::gmres ) ); ASSERT_EQ( "fgmres", toString( EnumType::fgmres ) ); ASSERT_EQ( "bicgstab", toString( EnumType::bicgstab ) ); + ASSERT_EQ( "richardson", toString( EnumType::richardson ) ); ASSERT_EQ( "preconditioner", toString( EnumType::preconditioner ) ); } @@ -48,9 +49,9 @@ TEST( LinearSolverParametersEnums, PreconditionerType ) ASSERT_EQ( "sgs", toString( EnumType::sgs ) ); ASSERT_EQ( "l1sgs", toString( EnumType::l1sgs ) ); ASSERT_EQ( "chebyshev", toString( EnumType::chebyshev ) ); - ASSERT_EQ( "iluk", toString( EnumType::iluk ) ); + ASSERT_EQ( "ilu", toString( EnumType::ilu ) ); ASSERT_EQ( "ilut", toString( EnumType::ilut ) ); - ASSERT_EQ( "icc", toString( EnumType::ic ) ); // Notice the discrepancy here + ASSERT_EQ( "ic", toString( EnumType::ic ) ); // Notice the discrepancy here ASSERT_EQ( "ict", toString( EnumType::ict ) ); ASSERT_EQ( "amg", toString( EnumType::amg ) ); ASSERT_EQ( "mgr", toString( EnumType::mgr ) ); @@ -133,9 +134,9 @@ TEST( LinearSolverParametersEnums, AMGSmootherType ) ASSERT_EQ( "sgs", toString( EnumType::sgs ) ); ASSERT_EQ( "l1sgs", toString( EnumType::l1sgs ) ); ASSERT_EQ( "chebyshev", toString( EnumType::chebyshev ) ); - ASSERT_EQ( "ilu0", toString( EnumType::ilu0 ) ); + ASSERT_EQ( "ilu", toString( EnumType::ilu ) ); ASSERT_EQ( "ilut", toString( EnumType::ilut ) ); - ASSERT_EQ( "ic0", toString( EnumType::ic0 ) ); + ASSERT_EQ( "ic", toString( EnumType::ic ) ); ASSERT_EQ( "ict", toString( EnumType::ict ) ); } diff --git a/src/coreComponents/linearAlgebra/utilities/BlockOperator.hpp b/src/coreComponents/linearAlgebra/utilities/BlockOperator.hpp index 3d91ea20501..732a1370414 100644 --- a/src/coreComponents/linearAlgebra/utilities/BlockOperator.hpp +++ b/src/coreComponents/linearAlgebra/utilities/BlockOperator.hpp @@ -57,17 +57,6 @@ class BlockOperator : public BlockOperatorView< VECTOR, OPERATOR > */ BlockOperator( BlockOperator const & rhs ); - /** - * @brief Move constructor. - * @param rhs the block operator to move from - */ - BlockOperator( BlockOperator && rhs ); - - /** - * @brief Destructor. - */ - virtual ~BlockOperator() override = default; - private: void setPointers(); @@ -106,14 +95,6 @@ BlockOperator< VECTOR, OPERATOR >::BlockOperator( BlockOperator const & rhs ) setPointers(); } -template< typename VECTOR, typename OPERATOR > -BlockOperator< VECTOR, OPERATOR >::BlockOperator( BlockOperator && rhs ) - : Base( std::move( rhs ) ), - m_operatorStorage( std::move( rhs.m_operatorStorage ) ) -{ - setPointers(); -} - } // namespace geos #endif //GEOS_LINEARALGEBRA_UTILITIES_BLOCKOPERATOR_HPP_ diff --git a/src/coreComponents/linearAlgebra/utilities/BlockOperatorView.hpp b/src/coreComponents/linearAlgebra/utilities/BlockOperatorView.hpp index 82b69dc0dbc..7d03c789e14 100644 --- a/src/coreComponents/linearAlgebra/utilities/BlockOperatorView.hpp +++ b/src/coreComponents/linearAlgebra/utilities/BlockOperatorView.hpp @@ -54,11 +54,6 @@ class BlockOperatorView : public LinearOperator< BlockVectorView< VECTOR > > */ ///@{ - /** - * @brief Destructor. - */ - virtual ~BlockOperatorView() override = default; - /** * @brief Deleted copy assignment. * @return not callable diff --git a/src/coreComponents/linearAlgebra/utilities/InverseNormalOperator.hpp b/src/coreComponents/linearAlgebra/utilities/InverseNormalOperator.hpp index 3ab698cefe8..fb722bdbd1f 100644 --- a/src/coreComponents/linearAlgebra/utilities/InverseNormalOperator.hpp +++ b/src/coreComponents/linearAlgebra/utilities/InverseNormalOperator.hpp @@ -52,7 +52,7 @@ class ExplicitTransposeInverse m_transposeSolver.apply( src, dst ); } -public: +private: Matrix m_transposeMatrix; Solver m_transposeSolver; diff --git a/src/coreComponents/linearAlgebra/utilities/LAIHelperFunctions.hpp b/src/coreComponents/linearAlgebra/utilities/LAIHelperFunctions.hpp index f92fa2ccfec..8299c2c460b 100644 --- a/src/coreComponents/linearAlgebra/utilities/LAIHelperFunctions.hpp +++ b/src/coreComponents/linearAlgebra/utilities/LAIHelperFunctions.hpp @@ -205,144 +205,67 @@ MATRIX permuteMatrix( MATRIX const & matrix, /** * @brief Computes rigid body modes * @tparam VECTOR output vector type - * @param mesh the mesh - * @param dofManager the degree-of-freedom manager - * @param selection list of field names - * @param rigidBodyModes the output array of linear algebra vectors containing RBMs + * @param nodePosition array of node coordinates + * @param dofIndex array of nodal degree-of-freedom indices + * @param dofOffset global dof offset for displacement field + * @param numLocalDof the number of locally owned displacement dofs + * @return the output array of linear algebra vectors containing RBMs */ template< typename VECTOR > -void computeRigidBodyModes( MeshLevel const & mesh, - DofManager const & dofManager, - std::vector< string > const & selection, - array1d< VECTOR > & rigidBodyModes ) +array1d< VECTOR > +computeRigidBodyModes( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, + arrayView1d< globalIndex const > const & dofIndex, + globalIndex const dofOffset, + localIndex const numLocalDof ) { - NodeManager const & nodeManager = mesh.getNodeManager(); + GEOS_ASSERT_EQ( nodePosition.size( 0 ), dofIndex.size() ); + integer const numComponents = nodePosition.size( 1 ); + integer const numRidigBodyModes = numComponents * ( numComponents + 1 ) / 2; - localIndex numComponents = 0; - array1d< localIndex > globalNodeList; - for( localIndex k = 0; k < LvArray::integerConversion< localIndex >( selection.size() ); ++k ) - { - if( dofManager.location( selection[k] ) == FieldLocation::Node ) - { - string const & dispDofKey = dofManager.getKey( selection[k] ); - arrayView1d< globalIndex const > const & dofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); - localIndex const numComponentsField = dofManager.numComponents( selection[k] ); - numComponents = numComponents > 0 ? numComponents : numComponentsField; - GEOS_ERROR_IF( numComponents != numComponentsField, "Rigid body modes called with different number of components." ); - globalIndex const globalOffset = dofManager.globalOffset( selection[k] ); - globalIndex const numLocalDofs = LvArray::integerConversion< globalIndex >( dofManager.numLocalDofs( selection[k] ) ); - for( globalIndex i = 0; i < dofNumber.size(); ++i ) - { - if( dofNumber[i] >= globalOffset && ( dofNumber[i] - globalOffset ) < numLocalDofs ) - { - globalNodeList.emplace_back( LvArray::integerConversion< localIndex >( dofNumber[i] - globalOffset ) / numComponentsField ); - } - } - } - } - localIndex const numNodes = globalNodeList.size(); - arrayView1d< localIndex const > globalNodeListView = globalNodeList.toViewConst(); - - arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition = nodeManager.referencePosition(); + array1d< VECTOR > rigidBodyModes( numRidigBodyModes ); - localIndex const numRidigBodyModes = numComponents * ( numComponents + 1 ) / 2; - rigidBodyModes.resize( numRidigBodyModes ); + // Translation RBMs for( localIndex k = 0; k < numComponents; ++k ) { - rigidBodyModes[k].create( numNodes * numComponents, MPI_COMM_GEOSX ); + rigidBodyModes[k].create( numLocalDof, MPI_COMM_GEOSX ); arrayView1d< real64 > const values = rigidBodyModes[k].open(); - forAll< parallelHostPolicy >( numNodes, [=]( localIndex const i ) + forAll< parallelHostPolicy >( dofIndex.size(), [=]( localIndex const i ) { - values[numComponents * i + k] = 1.0; + localIndex const localDof = LvArray::integerConversion< localIndex >( dofIndex[i] - dofOffset ); + if( 0 <= localDof && localDof < numLocalDof ) + { + values[localDof + k] = 1.0; + } } ); rigidBodyModes[k].close(); rigidBodyModes[k].scale( 1.0 / rigidBodyModes[k].norm2() ); } - switch( numComponents ) - { - case 2: - { - localIndex const k = 2; - rigidBodyModes[k].create( numNodes*numComponents, MPI_COMM_GEOSX ); - { - arrayView1d< real64 > const values = rigidBodyModes[k].open(); - forAll< parallelHostPolicy >( numNodes, [=]( localIndex const i ) - { - values[numComponents * i + 0] = -nodePosition[globalNodeListView[i]][1]; - values[numComponents * i + 1] = +nodePosition[globalNodeListView[i]][0]; - } ); - rigidBodyModes[k].close(); - } - for( localIndex j = 0; j < k; ++j ) - { - rigidBodyModes[k].axpy( -rigidBodyModes[k].dot( rigidBodyModes[j] ), rigidBodyModes[j] ); - } - rigidBodyModes[k].scale( 1.0 / rigidBodyModes[k].norm2() ); - break; - } - case 3: + // Rotation RBMs + for( localIndex k = numComponents; k < numRidigBodyModes; ++k ) + { + rigidBodyModes[k].create( numLocalDof, MPI_COMM_GEOSX ); + arrayView1d< real64 > const values = rigidBodyModes[k].open(); + integer const ind[2] = { ( k - numComponents + 1 ) % numComponents, + ( k - numComponents + 2 ) % numComponents }; + forAll< parallelHostPolicy >( dofIndex.size(), [=]( localIndex const i ) { - localIndex k = 3; - rigidBodyModes[k].create( numNodes*numComponents, MPI_COMM_GEOSX ); - { - arrayView1d< real64 > const values = rigidBodyModes[k].open(); - forAll< parallelHostPolicy >( numNodes, [=]( localIndex const i ) - { - values[numComponents * i + 0] = +nodePosition[globalNodeListView[i]][1]; - values[numComponents * i + 1] = -nodePosition[globalNodeListView[i]][0]; - } ); - rigidBodyModes[k].close(); - } - - for( localIndex j = 0; j < k; ++j ) - { - rigidBodyModes[k].axpy( -rigidBodyModes[k].dot( rigidBodyModes[j] ), rigidBodyModes[j] ); - } - rigidBodyModes[k].scale( 1.0 / rigidBodyModes[k].norm2() ); - - ++k; - rigidBodyModes[k].create( numNodes*numComponents, MPI_COMM_GEOSX ); - { - arrayView1d< real64 > const values = rigidBodyModes[k].open(); - forAll< parallelHostPolicy >( numNodes, [=]( localIndex const i ) - { - values[numComponents * i + 1] = -nodePosition[globalNodeListView[i]][2]; - values[numComponents * i + 2] = +nodePosition[globalNodeListView[i]][1]; - } ); - rigidBodyModes[k].close(); - } - - for( localIndex j = 0; j < k; ++j ) - { - rigidBodyModes[k].axpy( -rigidBodyModes[k].dot( rigidBodyModes[j] ), rigidBodyModes[j] ); - } - rigidBodyModes[k].scale( 1.0 / rigidBodyModes[k].norm2() ); - - ++k; - rigidBodyModes[k].create( numNodes*numComponents, MPI_COMM_GEOSX ); - { - arrayView1d< real64 > const values = rigidBodyModes[k].open(); - forAll< parallelHostPolicy >( numNodes, [=]( localIndex const i ) - { - values[numComponents * i + 0] = +nodePosition[globalNodeListView[i]][2]; - values[numComponents * i + 2] = -nodePosition[globalNodeListView[i]][0]; - } ); - rigidBodyModes[k].close(); - } - - for( localIndex j = 0; j < k; ++j ) + localIndex const localDof = LvArray::integerConversion< localIndex >( dofIndex[i] - dofOffset ); + if( 0 <= localDof && localDof < numLocalDof ) { - rigidBodyModes[k].axpy( -rigidBodyModes[k].dot( rigidBodyModes[j] ), rigidBodyModes[j] ); + values[localDof + ind[0]] = -nodePosition( i, ind[1] ); + values[localDof + ind[1]] = +nodePosition( i, ind[0] ); } - rigidBodyModes[k].scale( 1.0 / rigidBodyModes[k].norm2() ); - break; - } - default: + } ); + rigidBodyModes[k].close(); + for( localIndex j = 0; j < k; ++j ) { - GEOS_ERROR( "Rigid body modes computation unsupported for " << numComponents << " components." ); + rigidBodyModes[k].axpy( -rigidBodyModes[k].dot( rigidBodyModes[j] ), rigidBodyModes[j] ); } + rigidBodyModes[k].scale( 1.0 / rigidBodyModes[k].norm2() ); } + + return rigidBodyModes; } } // LAIHelperFunctions namespace diff --git a/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp b/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp index 46c6280262d..96efcc29070 100644 --- a/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp +++ b/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp @@ -20,6 +20,7 @@ #define GEOS_LINEARALGEBRA_UTILITIES_LINEARSOLVERPARAMETERS_HPP_ #include "codingUtilities/EnumStrings.hpp" +#include "common/DataTypes.hpp" namespace geos { @@ -42,6 +43,7 @@ struct LinearSolverParameters gmres, ///< GMRES fgmres, ///< Flexible GMRES bicgstab, ///< BiCGStab + richardson, ///< Richardson iteration preconditioner ///< Preconditioner only }; @@ -57,7 +59,7 @@ struct LinearSolverParameters sgs, ///< Symmetric Gauss-Seidel smoothing l1sgs, ///< l1-Symmetric Gauss-Seidel smoothing chebyshev, ///< Chebyshev polynomial smoothing - iluk, ///< Incomplete LU with k-level of fill + ilu, ///< Incomplete LU with k-level of fill ilut, ///< Incomplete LU with thresholding ic, ///< Incomplete Cholesky ict, ///< Incomplete Cholesky with thresholding @@ -65,7 +67,7 @@ struct LinearSolverParameters mgr, ///< Multigrid reduction (Hypre only) block, ///< Block preconditioner direct, ///< Direct solver as preconditioner - bgs, ///< Gauss-Seidel smoothing (backward sweep) + bgs ///< Gauss-Seidel smoothing (backward sweep) }; integer logLevel = 0; ///< Output level [0=none, 1=basic, 2=everything] @@ -74,7 +76,7 @@ struct LinearSolverParameters integer stopIfError = 1; ///< Whether to stop the simulation if the linear solver reports an error SolverType solverType = SolverType::direct; ///< Solver type - PreconditionerType preconditionerType = PreconditionerType::iluk; ///< Preconditioner type + PreconditionerType preconditionerType = PreconditionerType::ilu; ///< Preconditioner type /// Direct solver parameters: used for SuperLU_Dist interface through hypre and PETSc struct Direct @@ -122,6 +124,21 @@ struct LinearSolverParameters } krylov; ///< Krylov-method parameter struct + /// Relaxation/stationary iteration parameters (Richardson, damped Jacobi, etc.) + struct Relaxation + { + real64 weight = 2.0 / 3.0; ///< Relaxation weight (omega) for stationary iterations + } + relaxation; ///< Relaxation method parameters + + /// Chebyshev iteration/smoothing parameters + struct Chebyshev + { + integer order = 2; ///< Chebyshev order + integer eigNumIter = 10; ///< Number of eigenvalue estimation CG iterations + } + chebyshev; ///< Chebyshev smoother parameters + /// Matrix-scaling parameters struct Scaling { @@ -159,9 +176,9 @@ struct LinearSolverParameters sgs, ///< Symmetric Gauss-Seidel smoothing l1sgs, ///< l1-Symmetric Gauss-Seidel smoothing chebyshev, ///< Chebyshev polynomial smoothing - ilu0, ///< ILU(0) + ilu, ///< Incomplete LU ilut, ///< Incomplete LU with thresholding - ic0, ///< Incomplete Cholesky + ic, ///< Incomplete Cholesky ict ///< Incomplete Cholesky with thresholding }; @@ -236,6 +253,7 @@ struct LinearSolverParameters #endif integer maxLevels = 20; ///< Maximum number of coarsening levels + integer numCycles = 1; ///< Number of multigrid cycles CycleType cycleType = CycleType::V; ///< AMG cycle type CoarseType coarseType = CoarseType::direct; ///< Coarse-level solver/smoother InterpType interpolationType = InterpType::extendedI; ///< Interpolation algorithm @@ -245,6 +263,7 @@ struct LinearSolverParameters integer numFunctions = 1; ///< Number of amg functions integer aggressiveNumPaths = 1; ///< Number of paths agg. coarsening. integer aggressiveNumLevels = 0; ///< Number of levels for aggressive coarsening. + integer maxCoarseSize = 9; ///< Threshold for coarse grid size AggInterpType aggressiveInterpType = AggInterpType::multipass; ///< Interp. type for agg. coarsening. PreOrPost preOrPostSmoothing = PreOrPost::both; ///< Pre and/or post smoothing real64 threshold = 0.0; ///< Threshold for "strong connections" (for classical @@ -288,8 +307,7 @@ struct LinearSolverParameters StrategyType strategy = StrategyType::invalid; ///< Predefined MGR solution strategy (solver specific) integer separateComponents = false; ///< Apply a separate displacement component (SDC) filter before AMG construction string displacementFieldName; ///< Displacement field name need for SDC filter - integer areWellsShut = false; ///< Flag to let MGR know that wells are shut, and that jacobi can be applied to the - ///< well block + integer areWellsShut = false; ///< Flag to let MGR know that wells are shut and jacobi can be applied to the well block } mgr; ///< Multigrid reduction (MGR) parameters @@ -299,7 +317,7 @@ struct LinearSolverParameters integer fill = 0; ///< Fill level real64 threshold = 0.0; ///< Dropping threshold } - ifact; ///< Incomplete factorization parameter struct + ifact; ///< Incomplete factorization parameter struct /// Domain decomposition parameters struct DD @@ -307,6 +325,58 @@ struct LinearSolverParameters integer overlap = 0; ///< Ghost overlap } dd; ///< Domain decomposition parameter struct + + /// Block preconditioner parameters + struct Block + { + /// Shape of the block preconditioner + enum class Shape + { + Diagonal, ///< (D)^{-1} + UpperTriangular, ///< (DU)^{-1} + LowerTriangular, ///< (LD)^{-1} + LowerUpperTriangular ///< (LDU)^{-1} + }; + + /// Type of Schur complement approximation used + enum class SchurType + { + None, ///< No Schur complement - just block-GS/block-Jacobi preconditioner + FirstBlockDiagonal, ///< Approximate first block with its diagonal + RowsumDiagonalProbing, ///< Rowsum-preserving diagonal approximation constructed with probing + FirstBlockUserDefined ///< User defined preconditioner for the first block + }; + + /// Type of block row scaling to apply + enum class Scaling + { + None, ///< No scaling + FrobeniusNorm, ///< Equilibrate Frobenius norm of the diagonal blocks + UserProvided ///< User-provided scaling + }; + + Shape shape = Shape::UpperTriangular; ///< Block preconditioner shape + SchurType schurType = SchurType::RowsumDiagonalProbing; ///< Schur complement type + Scaling scaling = Scaling::FrobeniusNorm; ///< Type of system scaling to use + + array1d< LinearSolverParameters const * > subParams; ///< Pointers to parameters for sub-problems + array1d< integer > order; ///< Order of application of sub-problem solvers + + /** + * @brief Set the number of blocks and resize arrays accordingly. + * @param numBlocks the number of sub-problem blocks in the system + */ + void resize( integer const numBlocks ) + { + subParams.resize( numBlocks ); + order.resize( numBlocks ); + for( integer i = 0; i < numBlocks; ++i ) + { + order[i] = i; + } + } + } + block; ///< Block preconditioner parameters }; /// Declare strings associated with enumeration values. @@ -316,6 +386,7 @@ ENUM_STRINGS( LinearSolverParameters::SolverType, "gmres", "fgmres", "bicgstab", + "richardson", "preconditioner" ); /// Declare strings associated with enumeration values. @@ -327,9 +398,9 @@ ENUM_STRINGS( LinearSolverParameters::PreconditionerType, "sgs", "l1sgs", "chebyshev", - "iluk", + "ilu", "ilut", - "icc", + "ic", "ict", "amg", "mgr", @@ -396,9 +467,9 @@ ENUM_STRINGS( LinearSolverParameters::AMG::SmootherType, "sgs", "l1sgs", "chebyshev", - "ilu0", + "ilu", "ilut", - "ic0", + "ic", "ict" ); /// Declare strings associated with enumeration values. @@ -453,6 +524,26 @@ ENUM_STRINGS( LinearSolverParameters::AMG::NullSpaceType, "constantModes", "rigidBodyModes" ); +/// Declare strings associated with enumeration values. +ENUM_STRINGS( LinearSolverParameters::Block::Shape, + "D", + "DU", + "LD", + "LDU" ); + +/// Declare strings associated with enumeration values. +ENUM_STRINGS( LinearSolverParameters::Block::SchurType, + "none", + "diagonal", + "probing", + "user" ); + +/// Declare strings associated with enumeration values. +ENUM_STRINGS( LinearSolverParameters::Block::Scaling, + "none", + "frobenius", + "user" ); + } /* namespace geos */ #endif /*GEOS_LINEARALGEBRA_UTILITIES_LINEARSOLVERPARAMETERS_HPP_ */ diff --git a/src/coreComponents/linearAlgebra/utilities/NormalOperator.hpp b/src/coreComponents/linearAlgebra/utilities/NormalOperator.hpp index 39b22017893..01e9a2f5489 100644 --- a/src/coreComponents/linearAlgebra/utilities/NormalOperator.hpp +++ b/src/coreComponents/linearAlgebra/utilities/NormalOperator.hpp @@ -15,8 +15,8 @@ /** * @file NormalOperator.hpp */ -#ifndef GEOS_LINEARALGEBRA_NORMALOPERATOR_HPP_ -#define GEOS_LINEARALGEBRA_NORMALOPERATOR_HPP_ +#ifndef GEOS_LINEARALGEBRA_UTILITIES_NORMALOPERATOR_HPP_ +#define GEOS_LINEARALGEBRA_UTILITIES_NORMALOPERATOR_HPP_ #include "common/LinearOperator.hpp" @@ -27,19 +27,19 @@ namespace geos * @brief Wraps a matrix A and represents A^T * A as a linear operator. * @tparam LAI the linear algebra interface */ -template< typename LAI > -class NormalOperator : public LinearOperator< typename LAI::ParallelVector > +template< typename MATRIX > +class NormalOperator : public LinearOperator< typename MATRIX::Vector > { public: /// Alias for base type - using Base = LinearOperator< typename LAI::ParallelVector >; + using Base = LinearOperator< typename MATRIX::Vector >; /// Alias for vector type using Vector = typename Base::Vector; /// Alias for matrix type - using Matrix = typename LAI::ParallelMatrix; + using Matrix = MATRIX; /** * @brief Constructor @@ -49,11 +49,6 @@ class NormalOperator : public LinearOperator< typename LAI::ParallelVector > : m_matrix( mat ) {} - /** - * @brief Destructor. - */ - virtual ~NormalOperator() override = default; - /** * @brief Apply operator to a vector. * @param src input vector @@ -115,4 +110,4 @@ class NormalOperator : public LinearOperator< typename LAI::ParallelVector > } // namespace geos -#endif //GEOS_LINEARALGEBRA_NORMALOPERATOR_HPP_ +#endif //GEOS_LINEARALGEBRA_UTILITIES_NORMALOPERATOR_HPP_ diff --git a/src/coreComponents/linearAlgebra/utilities/TransposeOperator.hpp b/src/coreComponents/linearAlgebra/utilities/TransposeOperator.hpp index 8e5ce025f9b..36581f60a51 100644 --- a/src/coreComponents/linearAlgebra/utilities/TransposeOperator.hpp +++ b/src/coreComponents/linearAlgebra/utilities/TransposeOperator.hpp @@ -16,8 +16,8 @@ * @file TransposeOperator.hpp */ -#ifndef GEOS_LINEARALGEBRA_TRANSPOSEMATRIXOPERATOR_HPP_ -#define GEOS_LINEARALGEBRA_TRANSPOSEMATRIXOPERATOR_HPP_ +#ifndef GEOS_LINEARALGEBRA_UTILITIES_TRANSPOSEOPERATOR_HPP_ +#define GEOS_LINEARALGEBRA_UTILITIES_TRANSPOSEOPERATOR_HPP_ #include "linearAlgebra/common/LinearOperator.hpp" @@ -26,21 +26,21 @@ namespace geos /** * @brief Simple class that wraps a matrix and represents its transpose as a linear operator. - * @tparam LAI the linear algebra interface + * @tparam MATRIX the linear algebra matrix type */ -template< typename LAI > -class TransposeOperator : public LinearOperator< typename LAI::ParallelVector > +template< typename MATRIX > +class TransposeOperator : public LinearOperator< typename MATRIX::Vector > { public: /// Alias for base type - using Base = LinearOperator< typename LAI::ParallelVector >; + using Base = LinearOperator< typename MATRIX::Vector >; /// Alias for vector type using Vector = typename Base::Vector; /// Alias for matrix type - using Matrix = typename LAI::ParallelMatrix; + using Matrix = MATRIX; /** * @brief Constructor. @@ -51,11 +51,6 @@ class TransposeOperator : public LinearOperator< typename LAI::ParallelVector > m_matrix( mat ) { } - /** - * @brief Destructor. - */ - virtual ~TransposeOperator() override = default; - /** * @brief Apply operator to a vector. * @param src input vector (x) @@ -69,8 +64,7 @@ class TransposeOperator : public LinearOperator< typename LAI::ParallelVector > } /** - * @brief Get the number of global rows. - * @return Number of global rows in the operator. + * @brief @return the number of global rows. */ virtual globalIndex numGlobalRows() const override { @@ -78,8 +72,7 @@ class TransposeOperator : public LinearOperator< typename LAI::ParallelVector > } /** - * @brief Get the number of global columns. - * @return Number of global columns in the operator. + * @brief @return the number of global columns. */ virtual globalIndex numGlobalCols() const override { @@ -87,8 +80,7 @@ class TransposeOperator : public LinearOperator< typename LAI::ParallelVector > } /** - * @brief Get the number of local rows. - * @return Number of local rows in the operator. + * @brief @return the number of local rows. */ virtual localIndex numLocalRows() const override { @@ -96,8 +88,7 @@ class TransposeOperator : public LinearOperator< typename LAI::ParallelVector > } /** - * @brief Get the number of local columns. - * @return Number of local columns in the operator. + * @brief @return the number of local columns. */ virtual localIndex numLocalCols() const override { @@ -120,4 +111,4 @@ class TransposeOperator : public LinearOperator< typename LAI::ParallelVector > } -#endif //GEOS_LINEARALGEBRA_TRANSPOSEMATRIXOPERATOR_HPP_ +#endif //GEOS_LINEARALGEBRA_UTILITIES_TRANSPOSEOPERATOR_HPP_ diff --git a/src/coreComponents/mesh/CMakeLists.txt b/src/coreComponents/mesh/CMakeLists.txt index 93ae9104bdf..aa61d6e53cc 100644 --- a/src/coreComponents/mesh/CMakeLists.txt +++ b/src/coreComponents/mesh/CMakeLists.txt @@ -45,7 +45,6 @@ set( mesh_headers generators/InternalWellGenerator.hpp generators/InternalWellboreGenerator.hpp generators/MeshGeneratorBase.hpp - generators/ParMETISInterface.hpp generators/PartitionDescriptor.hpp generators/PrismUtilities.hpp generators/WellGeneratorBase.hpp @@ -109,7 +108,6 @@ set( mesh_sources generators/InternalWellGenerator.cpp generators/InternalWellboreGenerator.cpp generators/MeshGeneratorBase.cpp - generators/ParMETISInterface.cpp generators/WellGeneratorBase.cpp mpiCommunications/CommID.cpp mpiCommunications/CommunicationTools.cpp @@ -129,7 +127,7 @@ set( mesh_sources utilities/ComputationalGeometry.cpp ) -set( dependencyList ${parallelDeps} schema dataRepository constitutive finiteElement parmetis metis ) +set( dependencyList ${parallelDeps} schema dataRepository constitutive finiteElement ) if( ENABLE_VTK ) message(STATUS "Adding VTKMeshGenerator sources and headers") @@ -153,6 +151,12 @@ if( ENABLE_VTK ) endif() endif() +if( ENABLE_PARMETIS ) + set( mesh_headers ${mesh_headers} generators/ParMETISInterface.hpp ) + set( mesh_sources ${mesh_sources} generators/ParMETISInterface.cpp ) + set( dependencyList ${dependencyList} parmetis metis ) +endif() + if( ENABLE_SCOTCH ) set( mesh_headers ${mesh_headers} generators/PTScotchInterface.hpp ) set( mesh_sources ${mesh_sources} generators/PTScotchInterface.cpp ) diff --git a/src/coreComponents/mesh/DomainPartition.hpp b/src/coreComponents/mesh/DomainPartition.hpp index 7a41174df1c..cdc2ae10bff 100644 --- a/src/coreComponents/mesh/DomainPartition.hpp +++ b/src/coreComponents/mesh/DomainPartition.hpp @@ -267,6 +267,21 @@ class DomainPartition : public dataRepository::Group std::vector< NeighborCommunicator > const & getNeighbors() const { return m_neighbors; }; + /** + * @brief Get a list of neighbor ranks. + * @return Container of neighbor ranks. + */ + std::vector< int > getNeighborRanks() const + { + std::vector< int > ranks; + ranks.reserve( m_neighbors.size() ); + for( NeighborCommunicator const & neighbor : m_neighbors ) + { + ranks.push_back( neighbor.neighborRank() ); + } + return ranks; + } + private: /** diff --git a/src/coreComponents/mesh/ElementRegionManager.cpp b/src/coreComponents/mesh/ElementRegionManager.cpp index b967a2b9e76..5da2d682797 100644 --- a/src/coreComponents/mesh/ElementRegionManager.cpp +++ b/src/coreComponents/mesh/ElementRegionManager.cpp @@ -61,8 +61,7 @@ void ElementRegionManager::setMaxGlobalIndex() { m_localMaxGlobalIndex = std::max( m_localMaxGlobalIndex, subRegion.maxGlobalIndex() ); } ); - - m_maxGlobalIndex = MpiWrapper::max( m_localMaxGlobalIndex, MPI_COMM_GEOSX ); + ObjectManagerBase::setMaxGlobalIndex(); } diff --git a/src/coreComponents/physicsSolvers/LinearSolverParameters.cpp b/src/coreComponents/physicsSolvers/LinearSolverParameters.cpp index 6fa15ae76e4..99926dd5c03 100644 --- a/src/coreComponents/physicsSolvers/LinearSolverParameters.cpp +++ b/src/coreComponents/physicsSolvers/LinearSolverParameters.cpp @@ -20,8 +20,79 @@ namespace geos { + using namespace dataRepository; +/** + * @brief Register an input block handler for a struct. + * @tparam CHILD type of child Group that will handle input for @p BLOCK + * @tparam T type for which input handling is needed + * @param parent pointer to parent group + * @param key the group key (XML tag) for sub-block + * @param data the struct to handle input for + * @note blocks are assumed optional for now, but it may be changed, + * in which case InputFlags should be provided as a parameter. + */ +template< typename CHILD, typename T > +void registerInputBlock( Group * parent, char const * const key, T & data ) +{ + parent->registerGroup( key, std::make_unique< CHILD >( key, parent, data ) ).setInputFlags( InputFlags::OPTIONAL ); +} + +class BlockParametersInput final : public dataRepository::Group +{ +public: + + /// Constructor + BlockParametersInput( string const & name, + Group * const parent, + LinearSolverParameters::Block & params ); + + virtual Group * createChild( string const & childKey, string const & childName ) override + { + GEOS_UNUSED_VAR( childKey, childName ); + return nullptr; + } + + /// Keys appearing in XML + struct viewKeyStruct + { + static constexpr char const * shapeString() { return "shape"; } + static constexpr char const * schurTypeString() { return "schurType"; } + static constexpr char const * scalingString() { return "scaling"; } + }; + +private: + + LinearSolverParameters::Block & m_parameters; +}; + +BlockParametersInput::BlockParametersInput( string const & name, + Group * const parent, + LinearSolverParameters::Block & params ) + : + Group( name, parent ), + m_parameters( params ) +{ + registerWrapper( viewKeyStruct::shapeString(), &m_parameters.shape ). + setDefaultValue( m_parameters.shape ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Block preconditioner shape, options: " + "``" + EnumStrings< LinearSolverParameters::Block::Shape >::concat( "``, ``" ) + "``" ); + + registerWrapper( viewKeyStruct::schurTypeString(), &m_parameters.schurType ). + setDefaultValue( m_parameters.schurType ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Schur complement type, options: " + "``" + EnumStrings< LinearSolverParameters::Block::SchurType >::concat( "``, ``" ) + "``" ); + + registerWrapper( viewKeyStruct::scalingString(), &m_parameters.scaling ). + setDefaultValue( m_parameters.scaling ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Block scaling type, options: " + "``" + EnumStrings< LinearSolverParameters::Block::Scaling >::concat( "``, ``" ) + "``" ); +} + LinearSolverParametersInput::LinearSolverParametersInput( string const & name, Group * const parent ) : @@ -34,13 +105,13 @@ LinearSolverParametersInput::LinearSolverParametersInput( string const & name, setApplyDefaultValue( m_parameters.solverType ). setInputFlag( InputFlags::OPTIONAL ). setDescription( "Linear solver type. Available options are: " - "``" + EnumStrings< LinearSolverParameters::SolverType >::concat( "|" ) + "``" ); + "``" + EnumStrings< LinearSolverParameters::SolverType >::concat( "``, ``" ) + "``" ); registerWrapper( viewKeyStruct::preconditionerTypeString(), &m_parameters.preconditionerType ). setApplyDefaultValue( m_parameters.preconditionerType ). setInputFlag( InputFlags::OPTIONAL ). setDescription( "Preconditioner type. Available options are: " - "``" + EnumStrings< LinearSolverParameters::PreconditionerType >::concat( "|" ) + "``" ); + "``" + EnumStrings< LinearSolverParameters::PreconditionerType >::concat( "``, ``" ) + "``" ); registerWrapper( viewKeyStruct::stopIfErrorString(), &m_parameters.stopIfError ). setApplyDefaultValue( m_parameters.stopIfError ). @@ -61,13 +132,13 @@ LinearSolverParametersInput::LinearSolverParametersInput( string const & name, setApplyDefaultValue( m_parameters.direct.colPerm ). setInputFlag( InputFlags::OPTIONAL ). setDescription( "How to permute the columns. Available options are: " - "``" + EnumStrings< LinearSolverParameters::Direct::ColPerm >::concat( "|" ) + "``" ); + "``" + EnumStrings< LinearSolverParameters::Direct::ColPerm >::concat( "``, ``" ) + "``" ); registerWrapper( viewKeyStruct::directRowPermString(), &m_parameters.direct.rowPerm ). setApplyDefaultValue( m_parameters.direct.rowPerm ). setInputFlag( InputFlags::OPTIONAL ). setDescription( "How to permute the rows. Available options are: " - "``" + EnumStrings< LinearSolverParameters::Direct::RowPerm >::concat( "|" ) + "``" ); + "``" + EnumStrings< LinearSolverParameters::Direct::RowPerm >::concat( "``, ``" ) + "``" ); registerWrapper( viewKeyStruct::directReplTinyPivotString(), &m_parameters.direct.replaceTinyPivot ). setApplyDefaultValue( m_parameters.direct.replaceTinyPivot ). @@ -113,6 +184,26 @@ LinearSolverParametersInput::LinearSolverParametersInput( string const & name, setInputFlag( InputFlags::OPTIONAL ). setDescription( "Weakest-allowed tolerance for adaptive method" ); + registerWrapper( viewKeyStruct::relaxationWeightString(), &m_parameters.relaxation.weight ). + setApplyDefaultValue( m_parameters.relaxation.weight ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Relaxation weight (omega) for stationary iterations" ); + + registerWrapper( viewKeyStruct::chebyshevOrderString(), &m_parameters.chebyshev.order ). + setApplyDefaultValue( m_parameters.chebyshev.order ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Chebyshev order" ); + + registerWrapper( viewKeyStruct::chebyshevEigNumIterString(), &m_parameters.chebyshev.eigNumIter ). + setApplyDefaultValue( m_parameters.chebyshev.eigNumIter ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Number of eigenvalue estimation CG iterations" ); + + registerWrapper( viewKeyStruct::amgNumCyclesString(), &m_parameters.amg.numCycles ). + setApplyDefaultValue( m_parameters.amg.numCycles ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "AMG number of cycles" ); + registerWrapper( viewKeyStruct::amgNumSweepsString(), &m_parameters.amg.numSweeps ). setApplyDefaultValue( m_parameters.amg.numSweeps ). setInputFlag( InputFlags::OPTIONAL ). @@ -122,7 +213,7 @@ LinearSolverParametersInput::LinearSolverParametersInput( string const & name, setApplyDefaultValue( m_parameters.amg.smootherType ). setInputFlag( InputFlags::OPTIONAL ). setDescription( "AMG smoother type. Available options are: " - "``" + EnumStrings< LinearSolverParameters::AMG::SmootherType >::concat( "|" ) + "``" ); + "``" + EnumStrings< LinearSolverParameters::AMG::SmootherType >::concat( "``, ``" ) + "``" ); registerWrapper( viewKeyStruct::amgRelaxWeight(), &m_parameters.amg.relaxWeight ). setApplyDefaultValue( m_parameters.amg.relaxWeight ). @@ -133,7 +224,7 @@ LinearSolverParametersInput::LinearSolverParametersInput( string const & name, setApplyDefaultValue( m_parameters.amg.coarseType ). setInputFlag( InputFlags::OPTIONAL ). setDescription( "AMG coarsest level solver/smoother type. Available options are: " - "``" + EnumStrings< LinearSolverParameters::AMG::CoarseType >::concat( "|" ) + "``" ); + "``" + EnumStrings< LinearSolverParameters::AMG::CoarseType >::concat( "``, ``" ) + "``" ); registerWrapper( viewKeyStruct::amgCoarseningString(), &m_parameters.amg.coarseningType ). setApplyDefaultValue( m_parameters.amg.coarseningType ). @@ -173,6 +264,11 @@ LinearSolverParametersInput::LinearSolverParametersInput( string const & name, setDescription( "AMG aggressive interpolation algorithm. Available options are: " "``" + EnumStrings< LinearSolverParameters::AMG::AggInterpType >::concat( "|" ) + "``" ); + registerWrapper( viewKeyStruct::amgMaxCoarseSizeString(), &m_parameters.amg.maxCoarseSize ). + setApplyDefaultValue( m_parameters.amg.maxCoarseSize ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "AMG threshold for coarse grid size" ); + registerWrapper( viewKeyStruct::amgThresholdString(), &m_parameters.amg.threshold ). setApplyDefaultValue( m_parameters.amg.threshold ). setInputFlag( InputFlags::OPTIONAL ). @@ -186,8 +282,8 @@ LinearSolverParametersInput::LinearSolverParametersInput( string const & name, registerWrapper( viewKeyStruct::amgNullSpaceTypeString(), &m_parameters.amg.nullSpaceType ). setApplyDefaultValue( m_parameters.amg.nullSpaceType ). setInputFlag( InputFlags::OPTIONAL ). - setDescription( "AMG near null space approximation. Available options are:" - "``" + EnumStrings< LinearSolverParameters::AMG::NullSpaceType >::concat( "|" ) + "``" ); + setDescription( "AMG near null space approximation. Available options are: " + "``" + EnumStrings< LinearSolverParameters::AMG::NullSpaceType >::concat( "``, ``" ) + "``" ); registerWrapper( viewKeyStruct::iluFillString(), &m_parameters.ifact.fill ). setApplyDefaultValue( m_parameters.ifact.fill ). @@ -198,6 +294,8 @@ LinearSolverParametersInput::LinearSolverParametersInput( string const & name, setApplyDefaultValue( m_parameters.ifact.threshold ). setInputFlag( InputFlags::OPTIONAL ). setDescription( "ILU(T) threshold factor" ); + + registerInputBlock< BlockParametersInput >( this, groupKeyStruct::blockString(), m_parameters.block ); } void LinearSolverParametersInput::postProcessInput() @@ -259,6 +357,13 @@ void LinearSolverParametersInput::postProcessInput() // TODO input validation for other AMG parameters ? } +Group * LinearSolverParametersInput::createChild( string const & childKey, + string const & childName ) +{ + GEOS_UNUSED_VAR( childKey, childName ); + return nullptr; +} + REGISTER_CATALOG_ENTRY( Group, LinearSolverParametersInput, string const &, Group * const ) } // namespace geos diff --git a/src/coreComponents/physicsSolvers/LinearSolverParameters.hpp b/src/coreComponents/physicsSolvers/LinearSolverParameters.hpp index ced2aa2b62a..a3f80d68213 100644 --- a/src/coreComponents/physicsSolvers/LinearSolverParameters.hpp +++ b/src/coreComponents/physicsSolvers/LinearSolverParameters.hpp @@ -42,23 +42,17 @@ class LinearSolverParametersInput : public dataRepository::Group { public: - LinearSolverParametersInput() = delete; - /// Constructor LinearSolverParametersInput( string const & name, Group * const parent ); - /// Copy constructor - LinearSolverParametersInput( LinearSolverParametersInput && ) = default; - - /// Destructor - virtual ~LinearSolverParametersInput() override = default; - /// Catalog name static string catalogName() { return "LinearSolverParameters"; } /// Postprocessing of input virtual void postProcessInput() override; + virtual Group * createChild( string const & childKey, string const & childName ) override final; + LinearSolverParameters const & get() const { return m_parameters; } @@ -101,6 +95,16 @@ class LinearSolverParametersInput : public dataRepository::Group /// Krylov weakest tolerance key static constexpr char const * krylovWeakTolString() { return "krylovWeakestTol"; } + /// Relaxation weight key + static constexpr char const * relaxationWeightString() { return "relaxationWeight"; } + + /// Chebyshev order key + static constexpr char const * chebyshevOrderString() { return "chebyshevOrder"; } + /// Number of eigenvalue estimation CG iterations key + static constexpr char const * chebyshevEigNumIterString() { return "chebyshevEigNumIter"; } + + /// AMG number of sweeps key + static constexpr char const * amgNumCyclesString() { return "amgNumCycles"; } /// AMG number of sweeps key static constexpr char const * amgNumSweepsString() { return "amgNumSweeps"; } /// AMG smoother type key @@ -129,6 +133,8 @@ class LinearSolverParametersInput : public dataRepository::Group static constexpr char const * amgAggressiveInterpTypeString() { return "amgAggressiveInterpType"; } /// AMG separate components flag static constexpr char const * amgSeparateComponentsString() { return "amgSeparateComponents"; } + /// AMG coarse grid threshold key + static constexpr char const * amgMaxCoarseSizeString() { return "amgMaxCoarseSize"; } /// ILU fill key static constexpr char const * iluFillString() { return "iluFill"; } @@ -136,6 +142,12 @@ class LinearSolverParametersInput : public dataRepository::Group static constexpr char const * iluThresholdString() { return "iluThreshold"; } }; + /// Keys appearing in XML + struct groupKeyStruct + { + static constexpr char const * blockString() { return "Block"; } + }; + private: LinearSolverParameters m_parameters; diff --git a/src/coreComponents/physicsSolvers/SolverBase.cpp b/src/coreComponents/physicsSolvers/SolverBase.cpp index d3908f86b33..76e10eb61a9 100644 --- a/src/coreComponents/physicsSolvers/SolverBase.cpp +++ b/src/coreComponents/physicsSolvers/SolverBase.cpp @@ -355,7 +355,7 @@ real64 SolverBase::setNextDtBasedOnNewtonIter( real64 const & currentDt ) real64 SolverBase::linearImplicitStep( real64 const & time_n, real64 const & dt, - integer const GEOS_UNUSED_PARAM( cycleNumber ), + integer const cycleNumber, DomainPartition & domain ) { // call setup for physics solver. Pre step allocations etc. @@ -396,17 +396,11 @@ real64 SolverBase::linearImplicitStep( real64 const & time_n, m_assemblyCallback( m_localMatrix, std::move( localRhsCopy ) ); } - // TODO: Trilinos currently requires this, re-evaluate after moving to Tpetra-based solvers - if( m_precond ) - { - m_precond->clear(); - } - // Compose parallel LA matrix out of local matrix m_matrix.create( m_localMatrix.toViewConst(), m_dofManager.numLocalDofs(), MPI_COMM_GEOSX ); // Output the linear system matrix/rhs for debugging purposes - debugOutputSystem( 0.0, 0, 0, m_matrix, m_rhs ); + debugOutputSystem( time_n, cycleNumber, 0, m_matrix, m_rhs ); // Solve the linear system solveLinearSystem( m_dofManager, m_matrix, m_rhs, m_solution ); @@ -793,7 +787,8 @@ bool SolverBase::solveNonlinearSystem( real64 const & time_n, for( newtonIter = 0; newtonIter < maxNewtonIter; ++newtonIter ) { - GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " Attempt: {:2}, ConfigurationIter: {:2}, NewtonIter: {:2}", dtAttempt, configurationLoopIter, newtonIter ) ); + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " Attempt: {:2}, ConfigurationIter: {:2}, NewtonIter: {:2}", + dtAttempt, configurationLoopIter, newtonIter ) ); // zero out matrix/rhs before assembly m_localMatrix.zero(); @@ -832,14 +827,7 @@ bool SolverBase::solveNonlinearSystem( real64 const & time_n, // get residual norm real64 residualNorm = calculateResidualNorm( time_n, stepDt, domain, m_dofManager, m_rhs.values() ); - GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " ( R ) = ( {:4.2e} ) ; ", residualNorm ) ); - if( newtonIter > 0 ) - { - GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " Last LinSolve(iter,res) = ( {:3}, {:4.2e} ) ; ", - m_linearSolverResult.numIterations, - m_linearSolverResult.residualReduction ) ); - } // if the residual norm is less than the Newton tolerance we denote that we have // converged and break from the Newton loop immediately. @@ -917,12 +905,6 @@ bool SolverBase::solveNonlinearSystem( real64 const & time_n, krylovParams.relTolerance = eisenstatWalker( residualNorm, lastResidual, krylovParams.weakestTol ); } - // TODO: Trilinos currently requires this, re-evaluate after moving to Tpetra-based solvers - if( m_precond ) - { - m_precond->clear(); - } - // Compose parallel LA matrix/rhs out of local LA matrix/rhs // m_matrix.create( m_localMatrix.toViewConst(), m_dofManager.numLocalDofs(), MPI_COMM_GEOSX ); @@ -936,6 +918,10 @@ bool SolverBase::solveNonlinearSystem( real64 const & time_n, // Increment the solver statistics for reporting purposes m_solverStatistics.logNonlinearIteration( m_linearSolverResult.numIterations ); + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " Linear solve: ( iter, res ) = ( {:3}, {:4.2e} )", + m_linearSolverResult.numIterations, + m_linearSolverResult.residualReduction ) ); + // Output the linear system solution for debugging purposes debugOutputSolution( time_n, cycleNumber, newtonIter, m_solution ); @@ -1017,6 +1003,12 @@ void SolverBase::setupSystem( DomainPartition & domain, solution.create( dofManager.numLocalDofs(), MPI_COMM_GEOSX ); } +std::unique_ptr< PreconditionerBase< LAInterface > > +SolverBase::createPreconditioner( DomainPartition & GEOS_UNUSED_PARAM( domain ) ) const +{ + return LAInterface::createPreconditioner( m_linearSolverParameters.get() ); +} + void SolverBase::assembleSystem( real64 const GEOS_UNUSED_PARAM( time ), real64 const GEOS_UNUSED_PARAM( dt ), DomainPartition & GEOS_UNUSED_PARAM( domain ), @@ -1063,8 +1055,7 @@ void debugOutputLAObject( T const & obj, { if( toScreen ) { - string const frame( screenName.size() + 1, '=' ); - GEOS_LOG_RANK_0( frame << "\n" << screenName << ":\n" << frame ); + GEOS_LOG_RANK_0( GEOS_FMT( "{2:=>{1}}\n{0}:\n{2:=>{1}}", screenName, screenName.size() + 1, "" ) ); GEOS_LOG( obj ); } @@ -1072,7 +1063,7 @@ void debugOutputLAObject( T const & obj, { string const filename = GEOS_FMT( "{}_{:06}_{:02}.mtx", filePrefix.c_str(), cycleNumber, nonlinearIteration ); obj.write( filename, LAIOutputFormat::MATRIX_MARKET ); - GEOS_LOG_RANK_0( screenName << " written to " << filename ); + GEOS_LOG_RANK_0( GEOS_FMT( "{} written to {}", screenName, filename ) ); } } @@ -1151,7 +1142,11 @@ void SolverBase::solveLinearSystem( DofManager const & dofManager, } else { - m_precond->setup( matrix ); + { + Stopwatch timer; + m_precond->setup( matrix ); + GEOS_LOG_RANK_0_IF( params.logLevel >= 1, GEOS_FMT( "[Preconditioner] setup complete ({:.3f} s)", timer.elapsedTime() ) ); + } std::unique_ptr< KrylovSolver< ParallelVector > > solver = KrylovSolver< ParallelVector >::create( params, matrix, *m_precond ); solver->solve( rhs, solution ); m_linearSolverResult = solver->result(); diff --git a/src/coreComponents/physicsSolvers/SolverBase.hpp b/src/coreComponents/physicsSolvers/SolverBase.hpp index 02f8d812062..a231ba52b3b 100644 --- a/src/coreComponents/physicsSolvers/SolverBase.hpp +++ b/src/coreComponents/physicsSolvers/SolverBase.hpp @@ -310,6 +310,14 @@ class SolverBase : public ExecutableGroup ParallelVector & solution, bool const setSparsity = true ); + /** + * @brief Create a preconditioner for this solver's linear system. + * @param domain the domain containing the mesh and fields + * @return the newly created preconditioner object + */ + virtual std::unique_ptr< PreconditionerBase< LAInterface > > + createPreconditioner( DomainPartition & domain ) const; + /** * @brief function to assemble the linear system matrix and rhs * @param time the time at the beginning of the step diff --git a/src/coreComponents/physicsSolvers/contact/LagrangianContactSolver.cpp b/src/coreComponents/physicsSolvers/contact/LagrangianContactSolver.cpp index 057cc2f97aa..deb890e3d46 100644 --- a/src/coreComponents/physicsSolvers/contact/LagrangianContactSolver.cpp +++ b/src/coreComponents/physicsSolvers/contact/LagrangianContactSolver.cpp @@ -715,25 +715,28 @@ real64 LagrangianContactSolver::calculateResidualNorm( real64 const & GEOS_UNUSE return globalResidualNorm[2]; } -void LagrangianContactSolver::createPreconditioner( DomainPartition const & domain ) +std::unique_ptr< PreconditionerBase< LAInterface > > +LagrangianContactSolver::createPreconditioner( DomainPartition & domain ) const { - if( m_linearSolverParameters.get().preconditionerType == LinearSolverParameters::PreconditionerType::block ) + LinearSolverParameters const & linParams = m_linearSolverParameters.get(); + if( linParams.preconditionerType == LinearSolverParameters::PreconditionerType::block ) { // TODO: move among inputs (xml) string const leadingBlockApproximation = "blockJacobi"; LinearSolverParameters mechParams = m_solidSolver->getLinearSolverParameters(); - // Because of boundary conditions - mechParams.isSymmetric = false; std::unique_ptr< BlockPreconditioner< LAInterface > > precond; std::unique_ptr< PreconditionerBase< LAInterface > > tracPrecond; + LinearSolverParameters::Block blockParams; + blockParams.shape = LinearSolverParameters::Block::Shape::LowerUpperTriangular; + blockParams.scaling = LinearSolverParameters::Block::Scaling::UserProvided; + if( leadingBlockApproximation == "jacobi" ) { - precond = std::make_unique< BlockPreconditioner< LAInterface > >( BlockShapeOption::LowerUpperTriangular, - SchurComplementOption::FirstBlockDiagonal, - BlockScalingOption::UserProvided ); + blockParams.schurType = LinearSolverParameters::Block::SchurType::FirstBlockDiagonal; + precond = std::make_unique< BlockPreconditioner< LAInterface > >( blockParams ); // Using GEOSX implementation of Jacobi preconditioner // tracPrecond = std::make_unique< PreconditionerJacobi< LAInterface > >(); @@ -744,9 +747,8 @@ void LagrangianContactSolver::createPreconditioner( DomainPartition const & doma } else if( leadingBlockApproximation == "blockJacobi" ) { - precond = std::make_unique< BlockPreconditioner< LAInterface > >( BlockShapeOption::LowerUpperTriangular, - SchurComplementOption::FirstBlockUserDefined, - BlockScalingOption::UserProvided ); + blockParams.schurType = LinearSolverParameters::Block::SchurType::FirstBlockUserDefined; + precond = std::make_unique< BlockPreconditioner< LAInterface > >( blockParams ); tracPrecond = std::make_unique< PreconditionerBlockJacobi< LAInterface > >( mechParams.dofsPerNode ); } else @@ -759,31 +761,19 @@ void LagrangianContactSolver::createPreconditioner( DomainPartition const & doma { { contact::traction::key(), { 3, true } } }, std::move( tracPrecond ) ); - if( mechParams.amg.nullSpaceType == LinearSolverParameters::AMG::NullSpaceType::rigidBodyModes ) - { - if( m_solidSolver->getRigidBodyModes().empty() ) - { - MeshLevel const & mesh = domain.getMeshBody( 0 ).getBaseDiscretization(); - LAIHelperFunctions::computeRigidBodyModes( mesh, - m_dofManager, - { solidMechanics::totalDisplacement::key() }, - m_solidSolver->getRigidBodyModes() ); - } - } - // Preconditioner for the Schur complement: mechPrecond - std::unique_ptr< PreconditionerBase< LAInterface > > mechPrecond = LAInterface::createPreconditioner( mechParams, m_solidSolver->getRigidBodyModes() ); precond->setupBlock( 1, - { { solidMechanics::totalDisplacement::key(), { 3, true } } }, - std::move( mechPrecond ) ); + { { fields::solidMechanics::totalDisplacement::key(), { 3, true } } }, + m_solidSolver->createPreconditioner( domain ) ); - m_precond = std::move( precond ); + return precond; } else { - //TODO: Revisit this part such that is coherent across physics solver - //m_precond = LAInterface::createPreconditioner( m_linearSolverParameters.get() ); + // Unomment to use GEOSX's implementations of Krylov solvers instead of LA backend's + //return SolverBase::createPreconditioner( domain ); } + return {}; } void LagrangianContactSolver::computeRotationMatrices( DomainPartition & domain ) const diff --git a/src/coreComponents/physicsSolvers/contact/LagrangianContactSolver.hpp b/src/coreComponents/physicsSolvers/contact/LagrangianContactSolver.hpp index e9910a10ce0..dd72236ceb6 100644 --- a/src/coreComponents/physicsSolvers/contact/LagrangianContactSolver.hpp +++ b/src/coreComponents/physicsSolvers/contact/LagrangianContactSolver.hpp @@ -64,6 +64,9 @@ class LagrangianContactSolver : public ContactSolverBase ParallelVector & solution, bool const setSparsity = true ) override; + virtual std::unique_ptr< PreconditionerBase< LAInterface > > + createPreconditioner( DomainPartition & domain ) const override; + virtual void implicitStepSetup( real64 const & time_n, real64 const & dt, @@ -156,8 +159,6 @@ class LagrangianContactSolver : public ContactSolverBase real64 m_initialResidual[3] = {0.0, 0.0, 0.0}; - void createPreconditioner( DomainPartition const & domain ); - void computeFaceDisplacementJump( DomainPartition & domain ); virtual void setConstitutiveNames( ElementSubRegionBase & subRegion ) const override; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp index 10f6897529d..1e727242c8e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp @@ -354,12 +354,6 @@ void SinglePhaseBase::initializePostInitialConditionsPreSubGroups() MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { - FieldIdentifiers fieldsToBeSync; - fieldsToBeSync.addElementFields( { fields::flow::pressure::key() }, - regionNames ); - - CommunicationTools::getInstance().synchronizeFields( fieldsToBeSync, mesh, domain.getNeighbors(), false ); - // Moved the following part from ImplicitStepSetup to here since it only needs to be initialized once // They will be updated in applySystemSolution and ImplicitStepComplete, respectively mesh.getElemManager().forElementSubRegions< CellElementSubRegion, SurfaceElementSubRegion >( regionNames, [&]( localIndex const, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp index eecbd88e191..236910c18b8 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp @@ -23,13 +23,13 @@ #include "constitutive/permeability/PermeabilityFields.hpp" #include "constitutive/ConstitutivePassThru.hpp" #include "discretizationMethods/NumericalMethodsManager.hpp" -#include "mainInterface/ProblemManager.hpp" -#include "mesh/mpiCommunications/CommunicationTools.hpp" #include "finiteVolume/BoundaryStencil.hpp" #include "finiteVolume/FiniteVolumeManager.hpp" #include "finiteVolume/FluxApproximationBase.hpp" #include "fieldSpecification/FieldSpecificationManager.hpp" #include "fieldSpecification/AquiferBoundaryCondition.hpp" +#include "mainInterface/ProblemManager.hpp" +#include "mesh/mpiCommunications/CommunicationTools.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp" #include "physicsSolvers/fluidFlow/SinglePhaseBaseKernels.hpp" @@ -114,6 +114,25 @@ void SinglePhaseFVM< BASE >::setupSystem( DomainPartition & domain, solution, setSparsity ); + if( !m_precond && m_linearSolverParameters.get().solverType != LinearSolverParameters::SolverType::direct ) + { + m_precond = createPreconditioner( domain ); + } +} + +template< typename BASE > +std::unique_ptr< PreconditionerBase< LAInterface > > +SinglePhaseFVM< BASE >::createPreconditioner( DomainPartition & domain ) const +{ + LinearSolverParameters const & linParams = m_linearSolverParameters.get(); + GEOS_UNUSED_VAR( domain ); + switch( linParams.preconditionerType ) + { + default: + { + return LAInterface::createPreconditioner( linParams ); + } + } } template< typename BASE > diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp index d19c5942572..a9d494bdca8 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp @@ -57,6 +57,7 @@ class SinglePhaseFVM : public BASE using BASE::m_localMatrix; using BASE::m_linearSolverParameters; using BASE::m_nonlinearSolverParameters; + using BASE::m_precond; // Aliasing public/protected members/methods of FlowSolverBase so we don't // have to use this->member etc. @@ -131,6 +132,9 @@ class SinglePhaseFVM : public BASE ParallelVector & solution, bool const setSparsity = true ) override; + virtual std::unique_ptr< PreconditionerBase< LAInterface > > + createPreconditioner( DomainPartition & domain ) const override; + virtual void applyBoundaryConditions( real64 const time_n, real64 const dt, diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp index d3928354c6f..68cc620b7e0 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp @@ -42,7 +42,8 @@ using namespace fields; SinglePhasePoromechanics::SinglePhasePoromechanics( const string & name, Group * const parent ) : Base( name, parent ), - m_isThermal( 0 ) + m_isThermal( 0 ), + m_systemScaling( 0 ) { this->registerWrapper( viewKeyStruct::isThermalString(), &m_isThermal ). setApplyDefaultValue( 0 ). @@ -54,10 +55,16 @@ SinglePhasePoromechanics::SinglePhasePoromechanics( const string & name, setInputFlag( InputFlags::FALSE ). setDescription( "Flag to indicate that the solver is going to perform stress initialization" ); - m_linearSolverParameters.get().mgr.strategy = LinearSolverParameters::MGR::StrategyType::singlePhasePoromechanics; - m_linearSolverParameters.get().mgr.separateComponents = true; - m_linearSolverParameters.get().mgr.displacementFieldName = solidMechanics::totalDisplacement::key(); - m_linearSolverParameters.get().dofsPerNode = 3; + registerWrapper( viewKeyStruct::linearSystemScalingString(), &m_systemScaling ). + setInputFlag( InputFlags::OPTIONAL ). + setDefaultValue( m_systemScaling ). + setDescription( "Whether block system scaling should be performed" ); + + LinearSolverParameters & linParams = m_linearSolverParameters.get(); + linParams.mgr.strategy = LinearSolverParameters::MGR::StrategyType::singlePhasePoromechanics; + linParams.mgr.separateComponents = true; + linParams.mgr.displacementFieldName = solidMechanics::totalDisplacement::key(); + linParams.dofsPerNode = 3; } void SinglePhasePoromechanics::registerDataOnMesh( Group & meshBodies ) @@ -179,12 +186,13 @@ void SinglePhasePoromechanics::setupSystem( DomainPartition & domain, if( !m_precond && m_linearSolverParameters.get().solverType != LinearSolverParameters::SolverType::direct ) { - createPreconditioner(); + m_precond = createPreconditioner( domain ); } } void SinglePhasePoromechanics::initializePostInitialConditionsPreSubGroups() { + SolverBase::initializePostInitialConditionsPreSubGroups(); arrayView1d< string const > const & poromechanicsTargetRegionNames = @@ -207,17 +215,28 @@ void SinglePhasePoromechanics::initializePostInitialConditionsPreSubGroups() FlowSolverBase::viewKeyStruct::isThermalString(), flowSolver()->getName() ) ); isFlowThermal = m_isThermal; + using StrategyType = LinearSolverParameters::MGR::StrategyType; + LinearSolverParameters & linParams = m_linearSolverParameters.get(); if( m_isThermal ) { - m_linearSolverParameters.get().mgr.strategy = LinearSolverParameters::MGR::StrategyType::thermalSinglePhasePoromechanics; + linParams.mgr.strategy = StrategyType::thermalSinglePhasePoromechanics; } else { - if( flowSolver()->getLinearSolverParameters().mgr.strategy == LinearSolverParameters::MGR::StrategyType::singlePhaseHybridFVM ) + if( flowSolver()->getLinearSolverParameters().mgr.strategy == StrategyType::singlePhaseHybridFVM ) { - m_linearSolverParameters.get().mgr.strategy = LinearSolverParameters::MGR::StrategyType::hybridSinglePhasePoromechanics; + linParams.mgr.strategy = StrategyType::hybridSinglePhasePoromechanics; } } + + // Populate sub-block solver parameters for block preconditioner + linParams.block.resize( 2 ); + linParams.block.subParams[toUnderlying( SolverType::SolidMechanics )] = &solidMechanicsSolver()->getLinearSolverParameters(); + linParams.block.subParams[toUnderlying( SolverType::Flow )] = &flowSolver()->getLinearSolverParameters(); + + // For classical fixed-stress scheme the order must be mechanics-flow + linParams.block.order[toUnderlying( SolverType::SolidMechanics )] = 0; + linParams.block.order[toUnderlying( SolverType::Flow )] = 1; } void SinglePhasePoromechanics::assembleSystem( real64 const time_n, @@ -340,31 +359,79 @@ void SinglePhasePoromechanics::assembleElementBasedTerms( real64 const time_n, solidMechanicsSolver()->getMaxForce() = LvArray::math::max( mechanicsMaxForce, poromechanicsMaxForce ); } -void SinglePhasePoromechanics::createPreconditioner() +std::unique_ptr< PreconditionerBase< LAInterface > > +SinglePhasePoromechanics::createPreconditioner( DomainPartition & domain ) const { - if( m_linearSolverParameters.get().preconditionerType == LinearSolverParameters::PreconditionerType::block ) + LinearSolverParameters const & linParams = m_linearSolverParameters.get(); + switch( linParams.preconditionerType ) { - auto precond = std::make_unique< BlockPreconditioner< LAInterface > >( BlockShapeOption::UpperTriangular, - SchurComplementOption::RowsumDiagonalProbing, - BlockScalingOption::FrobeniusNorm ); - - auto mechPrecond = LAInterface::createPreconditioner( solidMechanicsSolver()->getLinearSolverParameters() ); - precond->setupBlock( 0, - { { solidMechanics::totalDisplacement::key(), { 3, true } } }, - std::make_unique< SeparateComponentPreconditioner< LAInterface > >( 3, std::move( mechPrecond ) ) ); + case LinearSolverParameters::PreconditionerType::block: + { + auto precond = std::make_unique< BlockPreconditioner< LAInterface > >( linParams.block ); - auto flowPrecond = LAInterface::createPreconditioner( flowSolver()->getLinearSolverParameters() ); - precond->setupBlock( 1, - { { flow::pressure::key(), { 1, true } } }, - std::move( flowPrecond ) ); + precond->setupBlock( linParams.block.order[toUnderlying( SolverType::SolidMechanics )], + { { solidMechanics::totalDisplacement::key(), { 3, true } } }, + solidMechanicsSolver()->createPreconditioner( domain ) ); + precond->setupBlock( linParams.block.order[toUnderlying( SolverType::Flow )], + { { SinglePhaseBase::viewKeyStruct::elemDofFieldString(), { 1, true } } }, + flowSolver()->createPreconditioner( domain ) ); - m_precond = std::move( precond ); + return precond; + } + default: + { + return SolverBase::createPreconditioner( domain ); + } } - else +} + +void SinglePhasePoromechanics::solveLinearSystem( DofManager const & dofManager, + ParallelMatrix & matrix, + ParallelVector & rhs, + ParallelVector & solution ) +{ + if( m_systemScaling ) { - //TODO: Revisit this part such that is coherent across physics solver - //m_precond = LAInterface::createPreconditioner( m_linearSolverParameters.get() ); + // Only compute this once and reuse for the entire simulation + if( !m_scalingVector.created() ) + { + // TODO: currently only handles displacement and cell pressure blocks, ignores face pressure in HybridFVM + DofManager::SubComponent fields[2]; + fields[toUnderlying( SolverType::SolidMechanics )] = { solidMechanics::totalDisplacement::key(), DofManager::CompMask{ 3, true } }; + fields[toUnderlying( SolverType::Flow )] = { SinglePhaseBase::viewKeyStruct::elemDofFieldString(), DofManager::CompMask{ 1, true } }; + + real64 norms[2]; + for( integer i = 0; i < 2; ++i ) + { + ParallelMatrix P, A; + dofManager.makeRestrictor( { fields[i] }, matrix.comm(), true, P ); + matrix.multiplyPtAP( P, A ); + norms[i] = A.normFrobenius(); + } + real64 const scale[2] = { std::min( norms[1] / norms[0], 1.0 ), std::min( norms[0] / norms[1], 1.0 ) }; + + m_scalingVector.create( rhs.localSize(), rhs.comm() ); + m_scalingVector.set( 1.0 ); + + localIndex offset = 0; + arrayView1d< real64 > const values = m_scalingVector.open(); + for( integer i = 0; i < 2; ++i ) + { + localIndex const numDof = dofManager.numLocalDofs( fields[i].fieldName ); + forAll< parallelDevicePolicy<> >( numDof, [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + values[offset + k] = scale[i]; + } ); + offset += numDof; + } + m_scalingVector.close(); + } + + matrix.leftScale( m_scalingVector ); + rhs.pointwiseProduct( m_scalingVector, rhs ); } + + SolverBase::solveLinearSystem( dofManager, matrix, rhs, solution ); } void SinglePhasePoromechanics::updateState( DomainPartition & domain ) diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp index 111a9f22719..bb2288119f0 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp @@ -55,9 +55,6 @@ class SinglePhasePoromechanics : public CoupledSolver< SinglePhaseBase, SinglePhasePoromechanics( const string & name, Group * const parent ); - /// Destructor for the class - ~SinglePhasePoromechanics() override {} - /** * @brief name of the node manager in the object catalog * @return string that contains the catalog name to generate a new SinglePhasePoromechanics object through the object catalog. @@ -108,6 +105,9 @@ class SinglePhasePoromechanics : public CoupledSolver< SinglePhaseBase, ParallelVector & solution, bool const setSparsity = true ) override; + virtual std::unique_ptr< PreconditionerBase< LAInterface > > + createPreconditioner( DomainPartition & domain ) const override; + virtual void assembleSystem( real64 const time, real64 const dt, DomainPartition & domain, @@ -115,6 +115,11 @@ class SinglePhasePoromechanics : public CoupledSolver< SinglePhaseBase, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) override; + virtual void solveLinearSystem( DofManager const & dofManager, + ParallelMatrix & matrix, + ParallelVector & rhs, + ParallelVector & solution ) override; + virtual void updateState( DomainPartition & domain ) override; /* @@ -138,6 +143,9 @@ class SinglePhasePoromechanics : public CoupledSolver< SinglePhaseBase, /// Flag to indicate that the solver is going to perform stress initialization constexpr static char const * performStressInitializationString() { return "performStressInitialization"; } + + /// Flag to enable block row scaling of the coupled system prior to the linear solve + constexpr static char const * linearSystemScalingString() { return "linearSystemScaling"; } }; protected: @@ -169,8 +177,6 @@ class SinglePhasePoromechanics : public CoupledSolver< SinglePhaseBase, */ void averageMeanStressIncrement( DomainPartition & domain ); - void createPreconditioner(); - template< typename CONSTITUTIVE_BASE, typename KERNEL_WRAPPER, typename ... PARAMS > @@ -185,6 +191,12 @@ class SinglePhasePoromechanics : public CoupledSolver< SinglePhaseBase, /// Flag to indicate that the solver is going to perform stress initialization integer m_performStressInitialization; + + // input flag indicating whether system scaling should be performed + integer m_systemScaling = 0; + + // A vector to perform row/block-wise linear system scaling + ParallelVector m_scalingVector; }; template< typename CONSTITUTIVE_BASE, diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index 0597ecd60bc..db218d70caf 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -26,14 +26,16 @@ #include "kernels/FixedStressThermoPoromechanics.hpp" #include "codingUtilities/Utilities.hpp" +#include "common/GEOS_RAJA_Interface.hpp" #include "constitutive/ConstitutiveManager.hpp" #include "constitutive/contact/ContactBase.hpp" -#include "common/GEOS_RAJA_Interface.hpp" #include "discretizationMethods/NumericalMethodsManager.hpp" #include "fieldSpecification/FieldSpecificationManager.hpp" #include "fieldSpecification/TractionBoundaryCondition.hpp" #include "finiteElement/FiniteElementDiscretizationManager.hpp" -#include "LvArray/src/output.hpp" +#include "finiteElement/FiniteElementDiscretizationManager.hpp" +#include "finiteElement/Kinematics.h" +#include "linearAlgebra/utilities/LAIHelperFunctions.hpp" #include "mainInterface/ProblemManager.hpp" #include "mesh/DomainPartition.hpp" #include "mesh/FaceElementSubRegion.hpp" @@ -115,24 +117,14 @@ SolidMechanicsLagrangianFEM::SolidMechanicsLagrangianFEM( const string & name, setInputFlag( InputFlags::FALSE ). setDescription( "The maximum force contribution in the problem domain." ); -} - -void SolidMechanicsLagrangianFEM::postProcessInput() -{ - SolverBase::postProcessInput(); - + // Set physics-dependent parameters for linear solver LinearSolverParameters & linParams = m_linearSolverParameters.get(); - linParams.isSymmetric = true; linParams.dofsPerNode = 3; - linParams.amg.separateComponents = true; -} -SolidMechanicsLagrangianFEM::~SolidMechanicsLagrangianFEM() -{ - // TODO Auto-generated destructor stub + // Must change default value to avoid being overwritten if attribute not specified in XML + m_linearSolverParameters.getWrapper< integer >( LinearSolverParametersInput::viewKeyStruct::amgSeparateComponentsString() ).setApplyDefaultValue( true ); } - void SolidMechanicsLagrangianFEM::registerDataOnMesh( Group & meshBodies ) { SolverBase::registerDataOnMesh( meshBodies ); @@ -999,7 +991,42 @@ void SolidMechanicsLagrangianFEM::setupSystem( DomainPartition & domain, sparsityPattern.compress(); localMatrix.assimilate< parallelDevicePolicy<> >( std::move( sparsityPattern ) ); + if( !m_precond && m_linearSolverParameters.get().solverType != LinearSolverParameters::SolverType::direct ) + { + m_precond = createPreconditioner( domain ); + } +} + +std::unique_ptr< PreconditionerBase< LAInterface > > +SolidMechanicsLagrangianFEM::createPreconditioner( DomainPartition & domain ) const +{ + LinearSolverParameters const & linParams = m_linearSolverParameters.get(); + if( linParams.amg.nullSpaceType == LinearSolverParameters::AMG::NullSpaceType::rigidBodyModes ) + { + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel const & mesh, auto const & ) + { + // The first target mesh body/level is used to compute RBMs + if( m_rigidBodyModes.empty() ) + { + NodeManager const & nodeManager = mesh.getNodeManager(); + arrayView1d< globalIndex const > const dofIndex = + nodeManager.getReference< array1d< globalIndex > >( m_dofManager.getKey( solidMechanics::totalDisplacement::key() ) ); + m_rigidBodyModes = LAIHelperFunctions::computeRigidBodyModes< ParallelVector >( nodeManager.referencePosition(), + dofIndex, + m_dofManager.rankOffset( solidMechanics::totalDisplacement::key() ), + m_dofManager.numLocalDofs( solidMechanics::totalDisplacement::key() ) ); + } + } ); + } + + switch( linParams.preconditionerType ) + { + default: + { + return LAInterface::createPreconditioner( linParams, m_rigidBodyModes ); + } + } } void SolidMechanicsLagrangianFEM::assembleSystem( real64 const GEOS_UNUSED_PARAM( time_n ), @@ -1185,7 +1212,7 @@ SolidMechanicsLagrangianFEM:: totalResidualNorm = std::max( residual, totalResidualNorm ); } ); - if( getLogLevel() >= 1 && logger::internal::rank==0 ) + if( getLogLevel() >= 1 && logger::internal::rank == 0 ) { std::cout << GEOS_FMT( " ( R{} ) = ( {:4.2e} ) ; ", coupledSolverAttributePrefix(), totalResidualNorm ); } diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp index f56799dd38e..af4466ad07c 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp @@ -65,18 +65,6 @@ class SolidMechanicsLagrangianFEM : public SolverBase SolidMechanicsLagrangianFEM( const string & name, Group * const parent ); - - SolidMechanicsLagrangianFEM( SolidMechanicsLagrangianFEM const & ) = delete; - SolidMechanicsLagrangianFEM( SolidMechanicsLagrangianFEM && ) = default; - - SolidMechanicsLagrangianFEM & operator=( SolidMechanicsLagrangianFEM const & ) = delete; - SolidMechanicsLagrangianFEM & operator=( SolidMechanicsLagrangianFEM && ) = delete; - - /** - * destructor - */ - virtual ~SolidMechanicsLagrangianFEM() override; - /** * @return The string that may be used to generate a new instance from the SolverBase::CatalogInterface::CatalogType */ @@ -124,6 +112,9 @@ class SolidMechanicsLagrangianFEM : public SolverBase ParallelVector & solution, bool const setSparsity = false ) override; + virtual std::unique_ptr< PreconditionerBase< LAInterface > > + createPreconditioner( DomainPartition & domain ) const override; + virtual void assembleSystem( real64 const time, real64 const dt, @@ -264,18 +255,7 @@ class SolidMechanicsLagrangianFEM : public SolverBase real64 & getMaxForce() { return m_maxForce; } - arrayView1d< ParallelVector > const & getRigidBodyModes() const - { - return m_rigidBodyModes; - } - - array1d< ParallelVector > & getRigidBodyModes() - { - return m_rigidBodyModes; - } - protected: - virtual void postProcessInput() override final; virtual void initializePostInitialConditionsPreSubGroups() override final; @@ -293,8 +273,8 @@ class SolidMechanicsLagrangianFEM : public SolverBase MPI_iCommData m_iComm; bool m_isFixedStressPoromechanicsUpdate; - /// Rigid body modes - array1d< ParallelVector > m_rigidBodyModes; + /// Rigid body modes; TODO remove mutable hack + mutable array1d< ParallelVector > m_rigidBodyModes; private: virtual void setConstitutiveNames( ElementSubRegionBase & subRegion ) const override; diff --git a/src/coreComponents/schema/docs/Block.rst b/src/coreComponents/schema/docs/Block.rst new file mode 100644 index 00000000000..38250e26b95 --- /dev/null +++ b/src/coreComponents/schema/docs/Block.rst @@ -0,0 +1,11 @@ + + +========= =========================================== ========= ============================================================================= +Name Type Default Description +========= =========================================== ========= ============================================================================= +scaling geos_LinearSolverParameters_Block_Scaling frobenius Block scaling type, options: ``none``, ``frobenius``, ``user`` +schurType geos_LinearSolverParameters_Block_SchurType probing Schur complement type, options: ``none``, ``diagonal``, ``probing``, ``user`` +shape geos_LinearSolverParameters_Block_Shape DU Block preconditioner shape, options: ``D``, ``DU``, ``LD``, ``LDU`` +========= =========================================== ========= ============================================================================= + + diff --git a/src/coreComponents/schema/docs/Block_other.rst b/src/coreComponents/schema/docs/Block_other.rst new file mode 100644 index 00000000000..adf1c1b8aec --- /dev/null +++ b/src/coreComponents/schema/docs/Block_other.rst @@ -0,0 +1,9 @@ + + +==== ==== ============================ +Name Type Description +==== ==== ============================ + (no documentation available) +==== ==== ============================ + + diff --git a/src/coreComponents/schema/docs/Hydrofracture.rst b/src/coreComponents/schema/docs/Hydrofracture.rst index 0165e3f988d..659df8d9c82 100644 --- a/src/coreComponents/schema/docs/Hydrofracture.rst +++ b/src/coreComponents/schema/docs/Hydrofracture.rst @@ -9,6 +9,7 @@ flowSolverName string required Name of the flow solver used by initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. isMatrixPoroelastic integer 0 (no description available) isThermal integer 0 Flag indicating whether the problem is thermal or not. Set isThermal="1" to enable the thermal coupling +linearSystemScaling integer 0 Whether block system scaling should be performed logLevel integer 0 Log level maxNumResolves integer 10 Value to indicate how many resolves may be executed to perform surface generation after the execution of flow and mechanics solver. name string required A name is required for any non-unique nodes diff --git a/src/coreComponents/schema/docs/LinearSolverParameters.rst b/src/coreComponents/schema/docs/LinearSolverParameters.rst index 67d62e850dd..8fbb3550f25 100644 --- a/src/coreComponents/schema/docs/LinearSolverParameters.rst +++ b/src/coreComponents/schema/docs/LinearSolverParameters.rst @@ -6,24 +6,28 @@ Name Type Def amgAggressiveCoarseningLevels integer 0 AMG number of levels for aggressive coarsening amgAggressiveCoarseningPaths integer 1 AMG number of paths for aggressive coarsening amgAggressiveInterpType geos_LinearSolverParameters_AMG_AggInterpType multipass AMG aggressive interpolation algorithm. Available options are: ``default\|extendedIStage2\|standardStage2\|extendedStage2\|multipass\|modifiedExtended\|modifiedExtendedI\|modifiedExtendedE\|modifiedMultipass`` -amgCoarseSolver geos_LinearSolverParameters_AMG_CoarseType direct AMG coarsest level solver/smoother type. Available options are: ``default\|jacobi\|l1jacobi\|fgs\|sgs\|l1sgs\|chebyshev\|direct\|bgs`` +amgCoarseSolver geos_LinearSolverParameters_AMG_CoarseType direct AMG coarsest level solver/smoother type. Available options are: ``default``, ``jacobi``, ``l1jacobi``, ``fgs``, ``sgs``, ``l1sgs``, ``chebyshev``, ``direct``, ``bgs`` amgCoarseningType geos_LinearSolverParameters_AMG_CoarseningType HMIS AMG coarsening algorithm. Available options are: ``default\|CLJP\|RugeStueben\|Falgout\|PMIS\|HMIS`` amgInterpolationMaxNonZeros integer 4 AMG interpolation maximum number of nonzeros per row amgInterpolationType geos_LinearSolverParameters_AMG_InterpType extendedI AMG interpolation algorithm. Available options are: ``default\|modifiedClassical\|direct\|multipass\|extendedI\|standard\|extended\|directBAMG\|modifiedExtended\|modifiedExtendedI\|modifiedExtendedE`` -amgNullSpaceType geos_LinearSolverParameters_AMG_NullSpaceType constantModes AMG near null space approximation. Available options are:``constantModes\|rigidBodyModes`` +amgMaxCoarseSize integer 9 AMG threshold for coarse grid size +amgNullSpaceType geos_LinearSolverParameters_AMG_NullSpaceType constantModes AMG near null space approximation. Available options are: ``constantModes``, ``rigidBodyModes`` +amgNumCycles integer 1 AMG number of cycles amgNumFunctions integer 1 AMG number of functions amgNumSweeps integer 1 AMG smoother sweeps amgRelaxWeight real64 1 AMG relaxation factor for the smoother amgSeparateComponents integer 0 AMG apply separate component filter for multi-variable problems -amgSmootherType geos_LinearSolverParameters_AMG_SmootherType l1sgs AMG smoother type. Available options are: ``default\|jacobi\|l1jacobi\|fgs\|bgs\|sgs\|l1sgs\|chebyshev\|ilu0\|ilut\|ic0\|ict`` +amgSmootherType geos_LinearSolverParameters_AMG_SmootherType l1sgs AMG smoother type. Available options are: ``default``, ``jacobi``, ``l1jacobi``, ``fgs``, ``bgs``, ``sgs``, ``l1sgs``, ``chebyshev``, ``ilu``, ``ilut``, ``ic``, ``ict`` amgThreshold real64 0 AMG strength-of-connection threshold +chebyshevEigNumIter integer 10 Number of eigenvalue estimation CG iterations +chebyshevOrder integer 2 Chebyshev order directCheckResidual integer 0 Whether to check the linear system solution residual -directColPerm geos_LinearSolverParameters_Direct_ColPerm metis How to permute the columns. Available options are: ``none\|MMD_AtplusA\|MMD_AtA\|colAMD\|metis\|parmetis`` +directColPerm geos_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 geos_LinearSolverParameters_Direct_RowPerm mc64 How to permute the rows. Available options are: ``none\|mc64`` +directRowPerm geos_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 @@ -35,9 +39,11 @@ krylovTol real64 1e- | :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 geos_LinearSolverParameters_PreconditionerType iluk Preconditioner type. Available options are: ``none\|jacobi\|l1jacobi\|fgs\|sgs\|l1sgs\|chebyshev\|iluk\|ilut\|icc\|ict\|amg\|mgr\|block\|direct\|bgs`` -solverType geos_LinearSolverParameters_SolverType direct Linear solver type. Available options are: ``direct\|cg\|gmres\|fgmres\|bicgstab\|preconditioner`` +preconditionerType geos_LinearSolverParameters_PreconditionerType ilu Preconditioner type. Available options are: ``none``, ``jacobi``, ``l1jacobi``, ``fgs``, ``sgs``, ``l1sgs``, ``chebyshev``, ``ilu``, ``ilut``, ``ic``, ``ict``, ``amg``, ``mgr``, ``block``, ``direct``, ``bgs`` +relaxationWeight real64 0.666667 Relaxation weight (omega) for stationary iterations +solverType geos_LinearSolverParameters_SolverType direct Linear solver type. Available options are: ``direct``, ``cg``, ``gmres``, ``fgmres``, ``bicgstab``, ``richardson``, ``preconditioner`` stopIfError integer 1 Whether to stop the simulation if the linear solver reports an error +Block node unique :ref:`XML_Block` ============================= ============================================== ============= ======================================================================================================================================================================================================================================================================================================================= diff --git a/src/coreComponents/schema/docs/LinearSolverParameters_other.rst b/src/coreComponents/schema/docs/LinearSolverParameters_other.rst index adf1c1b8aec..ac9d234d064 100644 --- a/src/coreComponents/schema/docs/LinearSolverParameters_other.rst +++ b/src/coreComponents/schema/docs/LinearSolverParameters_other.rst @@ -1,9 +1,9 @@ -==== ==== ============================ -Name Type Description -==== ==== ============================ - (no documentation available) -==== ==== ============================ +===== ==== ========================== +Name Type Description +===== ==== ========================== +Block node :ref:`DATASTRUCTURE_Block` +===== ==== ========================== diff --git a/src/coreComponents/schema/docs/SinglePhasePoromechanics.rst b/src/coreComponents/schema/docs/SinglePhasePoromechanics.rst index 284c642ef53..21b488f573e 100644 --- a/src/coreComponents/schema/docs/SinglePhasePoromechanics.rst +++ b/src/coreComponents/schema/docs/SinglePhasePoromechanics.rst @@ -7,6 +7,7 @@ cflFactor real64 0.5 Factor to apply to the `CFL cond flowSolverName string required Name of the flow solver used by the coupled solver initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. isThermal integer 0 Flag indicating whether the problem is thermal or not. Set isThermal="1" to enable the thermal coupling +linearSystemScaling integer 0 Whether block system scaling should be performed logLevel integer 0 Log level name string required A name is required for any non-unique nodes solidSolverName string required Name of the solid solver used by the coupled solver diff --git a/src/coreComponents/schema/docs/SinglePhasePoromechanicsEmbeddedFractures.rst b/src/coreComponents/schema/docs/SinglePhasePoromechanicsEmbeddedFractures.rst index 1f0093c2686..8b5aab05242 100644 --- a/src/coreComponents/schema/docs/SinglePhasePoromechanicsEmbeddedFractures.rst +++ b/src/coreComponents/schema/docs/SinglePhasePoromechanicsEmbeddedFractures.rst @@ -8,6 +8,7 @@ flowSolverName string required Name of the flow solver used by fracturesSolverName string required Name of the fractures solver to use in the fractured poroelastic solver initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. isThermal integer 0 Flag indicating whether the problem is thermal or not. Set isThermal="1" to enable the thermal coupling +linearSystemScaling integer 0 Whether block system scaling should be performed logLevel integer 0 Log level name string required A name is required for any non-unique nodes solidSolverName string required Name of the solid solver used by the coupled solver diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index 6e073e4a00d..f6b4ed49a19 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -1555,13 +1555,16 @@ stress - traction is applied to the faces as specified by the inner product of i + + + - + @@ -1569,8 +1572,12 @@ stress - traction is applied to the faces as specified by the inner product of i - + + + + + @@ -1579,13 +1586,17 @@ stress - traction is applied to the faces as specified by the inner product of i - + + + + + - + @@ -1595,7 +1606,7 @@ stress - traction is applied to the faces as specified by the inner product of i - + @@ -1616,13 +1627,38 @@ the relative residual norm satisfies: - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -1650,7 +1686,7 @@ the relative residual norm satisfies: - + @@ -1665,12 +1701,12 @@ the relative residual norm satisfies: - + - + @@ -2410,6 +2446,8 @@ Equal to 1 for surface conditions, and to 0 for reservoir conditions--> + + @@ -2713,6 +2751,8 @@ Local - Add stabilization only to interiors of macro elements.--> + + @@ -2759,6 +2799,8 @@ Local - Add stabilization only to interiors of macro elements.--> + + diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other index d6e5a858135..8bc6860de4c 100644 --- a/src/coreComponents/schema/schema.xsd.other +++ b/src/coreComponents/schema/schema.xsd.other @@ -406,7 +406,12 @@ - + + + + + + diff --git a/src/docs/doxygen/GeosxConfig.hpp b/src/docs/doxygen/GeosxConfig.hpp index 98135d0b29e..6dcce547ead 100644 --- a/src/docs/doxygen/GeosxConfig.hpp +++ b/src/docs/doxygen/GeosxConfig.hpp @@ -89,6 +89,12 @@ /// Enables use of PETSc library (CMake option ENABLE_PETSC) #define GEOSX_USE_PETSC +/// Enables use of METIS library (CMake option ENABLE_METIS) +#define GEOSX_USE_METIS + +/// Enables use of ParMETIS library (CMake option ENABLE_PARMETIS) +#define GEOSX_USE_PARMETIS + /// Enables use of Scotch library (CMake option ENABLE_SCOTCH) #define GEOSX_USE_SCOTCH diff --git a/src/docs/sphinx/CompleteXMLSchema.rst b/src/docs/sphinx/CompleteXMLSchema.rst index eb52d4dbff6..a9326eca32d 100644 --- a/src/docs/sphinx/CompleteXMLSchema.rst +++ b/src/docs/sphinx/CompleteXMLSchema.rst @@ -50,6 +50,13 @@ Element: BlackOilFluid .. include:: ../../coreComponents/schema/docs/BlackOilFluid.rst +.. _XML_Block: + +Element: Block +============== +.. include:: ../../coreComponents/schema/docs/Block.rst + + .. _XML_Blueprint: Element: Blueprint @@ -1279,6 +1286,13 @@ Datastructure: BlackOilFluid .. include:: ../../coreComponents/schema/docs/BlackOilFluid_other.rst +.. _DATASTRUCTURE_Block: + +Datastructure: Block +==================== +.. include:: ../../coreComponents/schema/docs/Block_other.rst + + .. _DATASTRUCTURE_Blueprint: Datastructure: Blueprint