diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json index 9853a691e0c..d8ebb7ebf87 100644 --- a/.devcontainer/devcontainer.json +++ b/.devcontainer/devcontainer.json @@ -2,7 +2,7 @@ "build": { "dockerfile": "Dockerfile", "args": { - "GEOS_TPL_TAG": "312-717" + "GEOS_TPL_TAG": "314-735" } }, "runArgs": [ diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 70459eee6ee..be3911f4273 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3682-12632-86ad358 + baseline: integratedTests/baseline_integratedTests-pr2427-12699-941dbff allow_fail: all: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index f04892b1613..d14911130c7 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines. Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining. These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD). +PR #2427 (2025-08-09) +===================== +Change default value of amgNumFunctions from 1 to 3 for solid mechanics solvers. Change in mesh partitioning of PoroElastic_hybridHexPrism_co2 cases due to Scotch version update. + PR #3682 (2025-08-08) ===================== Add physics-based scaling option. diff --git a/inputFiles/poromechanics/PoroElastic_Mandel_smoke_sequential.xml b/inputFiles/poromechanics/PoroElastic_Mandel_smoke_sequential.xml index badeed3ad15..187f5256212 100644 --- a/inputFiles/poromechanics/PoroElastic_Mandel_smoke_sequential.xml +++ b/inputFiles/poromechanics/PoroElastic_Mandel_smoke_sequential.xml @@ -17,9 +17,9 @@ scotch_VERSION = ${SCOTCH_VERSION}") + + foreach(_scotch_target SCOTCH::scotch SCOTCH::ptscotch SCOTCH::scotcherr ) + get_target_property(SCOTCH_INCLUDE_DIRS ${_scotch_target} INTERFACE_INCLUDE_DIRECTORIES) + set_target_properties(${_scotch_target} PROPERTIES INTERFACE_SYSTEM_INCLUDE_DIRECTORIES "${SCOTCH_INCLUDE_DIRS}") + set(thirdPartyLibs ${thirdPartyLibs} ${_scotch_target}) + endforeach() else() if(ENABLE_SCOTCH) message(WARNING "ENABLE_SCOTCH is ON but SCOTCH_DIR isn't defined.") diff --git a/src/coreComponents/codingUtilities/UnitTestUtilities.hpp b/src/coreComponents/codingUtilities/UnitTestUtilities.hpp index 092a524e6c0..409a838f851 100644 --- a/src/coreComponents/codingUtilities/UnitTestUtilities.hpp +++ b/src/coreComponents/codingUtilities/UnitTestUtilities.hpp @@ -122,13 +122,11 @@ inline void checkRelativeError( real64 const v1, real64 const v2, real64 const r 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 ) @@ -153,17 +151,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 > @@ -178,45 +190,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 e971431e6da..8d7766fa914 100644 --- a/src/coreComponents/common/DataTypes.hpp +++ b/src/coreComponents/common/DataTypes.hpp @@ -301,12 +301,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 c015f3e8856..365749d1ea4 100644 --- a/src/coreComponents/common/GeosxConfig.hpp.in +++ b/src/coreComponents/common/GeosxConfig.hpp.in @@ -98,6 +98,12 @@ /// Enables use of PETSc library (CMake option ENABLE_PETSC) #cmakedefine GEOS_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 GEOS_USE_SCOTCH diff --git a/src/coreComponents/common/MpiWrapper.hpp b/src/coreComponents/common/MpiWrapper.hpp index b884e15a8bc..95613986dd4 100644 --- a/src/coreComponents/common/MpiWrapper.hpp +++ b/src/coreComponents/common/MpiWrapper.hpp @@ -382,6 +382,11 @@ struct MpiWrapper array1d< T > & recvbuf, MPI_Comm comm = MPI_COMM_GEOS ); + template< typename T > + static int allGatherv( arrayView1d< T const > const & sendbuf, + array1d< T > & recvbuf, + MPI_Comm comm = MPI_COMM_GEOS ); + /** * @brief Convenience wrapper for the MPI_Allreduce function. * @tparam T type of data to reduce. Must correspond to a valid MPI_Datatype. @@ -619,6 +624,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. @@ -985,6 +1006,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 GEOS_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, @@ -1239,6 +1292,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 GEOS_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 3d42f1d83b4..a948f6555b3 100644 --- a/src/coreComponents/common/TimingMacros.hpp +++ b/src/coreComponents/common/TimingMacros.hpp @@ -35,7 +35,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 ); } } @@ -96,10 +96,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 1837effe1f8..605f5395222 100644 --- a/src/coreComponents/linearAlgebra/CMakeLists.txt +++ b/src/coreComponents/linearAlgebra/CMakeLists.txt @@ -44,6 +44,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 @@ -70,6 +72,7 @@ set( linearAlgebra_sources solvers/CgSolver.cpp solvers/GmresSolver.cpp solvers/KrylovSolver.cpp + solvers/RichardsonSolver.cpp solvers/SeparateComponentPreconditioner.cpp utilities/ReverseCutHillMcKeeOrdering.cpp ) diff --git a/src/coreComponents/linearAlgebra/DofManager.cpp b/src/coreComponents/linearAlgebra/DofManager.cpp index 4625bca014c..0a329859c4c 100644 --- a/src/coreComponents/linearAlgebra/DofManager.cpp +++ b/src/coreComponents/linearAlgebra/DofManager.cpp @@ -67,11 +67,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 ); } @@ -149,6 +149,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 { @@ -334,7 +339,7 @@ void DofManager::removeIndexArray( FieldDescription const & field ) } ); } -void DofManager::computeFieldDimensions( localIndex const fieldIndex ) +void DofManager::computeFieldDimensions( integer const fieldIndex ) { FieldDescription & field = m_fields[fieldIndex]; @@ -436,7 +441,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" ); stdVector< FieldSupport > processedSupports = processFieldSupportList( *m_domain, regions ); @@ -606,8 +611,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 stdVector< FieldSupport > processSupportList = processCouplingRegionList( supports, @@ -650,7 +655,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" ); @@ -830,11 +835,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; @@ -876,7 +881,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; } @@ -905,11 +910,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() ); } @@ -918,8 +923,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 ); @@ -1045,14 +1050,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; @@ -1100,7 +1105,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; } @@ -1111,8 +1116,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 ); @@ -1240,13 +1245,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 ); } @@ -1256,9 +1261,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 ); } @@ -1582,7 +1587,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() ); } @@ -1596,7 +1601,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; @@ -1665,6 +1670,7 @@ void DofManager::setupFrom( DofManager const & source, stdVector< 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(), @@ -1682,8 +1688,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 } ); } } @@ -1699,7 +1705,6 @@ void DofManager::makeRestrictor( stdVector< 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() ); stdVector< FieldDescription > fieldsSelected( selection.size() ); for( std::size_t k = 0; k < fieldsSelected.size(); ++k ) @@ -1724,7 +1729,7 @@ void DofManager::makeRestrictor( stdVector< 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; } ); @@ -1736,7 +1741,7 @@ void DofManager::makeRestrictor( stdVector< 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; } ); @@ -1748,44 +1753,44 @@ void DofManager::makeRestrictor( stdVector< 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(); } @@ -1795,12 +1800,12 @@ void DofManager::printFieldInfo( std::ostream & os ) const { if( MpiWrapper::commRank( MPI_COMM_GEOS ) == 0 ) { - localIndex const numFields = LvArray::integerConversion< localIndex >( m_fields.size() ); + integer const numFields = LvArray::integerConversion< integer >( m_fields.size() ); TableLayout const fieldLayout( "The summary of declared fields and coupling", {" ", "name", "comp", "N global DOF"} ); TableData fieldData; - for( localIndex i = 0; i < numFields; ++i ) + for( integer i = 0; i < numFields; ++i ) { FieldDescription const & f = m_fields[i]; fieldData.addRow( i, f.name, f.numComponents, f.numGlobalDof ); @@ -1809,9 +1814,9 @@ void DofManager::printFieldInfo( std::ostream & os ) const GEOS_LOG_RANK_0( logFormatter.toString( fieldData ) ); 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 ) { @@ -1862,7 +1867,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 87c36f1d1d0..73111b8727a 100644 --- a/src/coreComponents/linearAlgebra/DofManager.hpp +++ b/src/coreComponents/linearAlgebra/DofManager.hpp @@ -21,6 +21,7 @@ #define GEOS_LINEARALGEBRA_DOFMANAGER_HPP_ #include "common/DataTypes.hpp" +#include "common/Span.hpp" #include "linearAlgebra/utilities/ComponentMask.hpp" #include "mesh/FieldIdentifiers.hpp" @@ -45,10 +46,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. @@ -346,6 +347,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 @@ -368,7 +375,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 ) { @@ -397,7 +404,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. @@ -413,7 +420,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. @@ -428,7 +435,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. @@ -443,7 +450,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 @@ -527,13 +534,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 @@ -573,11 +580,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. @@ -588,15 +595,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 @@ -635,14 +638,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 stdVector< 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 57476c305bc..f3e6eb962ee 100644 --- a/src/coreComponents/linearAlgebra/common/LinearOperator.hpp +++ b/src/coreComponents/linearAlgebra/common/LinearOperator.hpp @@ -39,12 +39,12 @@ class LinearOperator using Vector = VECTOR; /** - * @brief Constructor + * @brief Constructor. */ LinearOperator() = default; /** - * @brief Destructor + * @brief Destructor. */ virtual ~LinearOperator() = default; @@ -74,26 +74,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), @@ -101,6 +111,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 ac9d0b022a9..5a7cea8d7f8 100644 --- a/src/coreComponents/linearAlgebra/common/PreconditionerBase.hpp +++ b/src/coreComponents/linearAlgebra/common/PreconditionerBase.hpp @@ -35,8 +35,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 42a1b8db354..237c25bc9c3 100644 --- a/src/coreComponents/linearAlgebra/interfaces/MatrixBase.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/MatrixBase.hpp @@ -46,7 +46,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 @@ -85,10 +85,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 @@ -617,6 +619,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 multiplyP1tAP2( 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 @@ -625,14 +648,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 ); + multiplyP1tAP2( P, P, dst ); } /** @@ -784,13 +800,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. @@ -806,6 +833,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 @@ -822,6 +856,55 @@ 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() ); + rowLengths.move( LvArray::MemorySpace::host ); + 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 @@ -863,18 +946,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 53a45e2ce7f..33d3affe7cf 100644 --- a/src/coreComponents/linearAlgebra/interfaces/direct/SuiteSparse.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/direct/SuiteSparse.cpp @@ -21,6 +21,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" @@ -138,7 +139,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 ); } @@ -160,7 +161,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 ); } @@ -170,7 +171,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; } @@ -217,6 +218,7 @@ template< typename LAI > void SuiteSparse< LAI >::apply( Vector const & src, Vector & dst ) const { + GEOS_MARK_FUNCTION; doSolve( src, dst, false ); } @@ -272,7 +274,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() ) == 0 ) { GEOS_WARNING( "SuiteSparse: failed to reduce residual below tolerance.\n" "Condition number estimate: " << condEst ); @@ -282,14 +284,9 @@ void SuiteSparse< LAI >::solve( Vector const & rhs, } } - if( m_params.logLevel >= 1 ) - { - GEOS_LOG_RANK_0( " Linear Solver | " << m_result.status << - " | Iterations: " << m_result.numIterations << - " | Final Rel Res: " << m_result.residualReduction << - " | Setup Time: " << m_result.setupTime << " s" << - " | Solve Time: " << m_result.solveTime << " s" ); - } + GEOS_LOG_RANK_0_IF( m_params.logLevel >= 1, + GEOS_FMT( " Linear Solver | {} | Iterations: {} | Final Rel Res: {} | Setup Time: {} s | Solve Time: {} s", + m_result.status, m_result.numIterations, m_result.residualReduction, m_result.setupTime, m_result.solveTime ) ); } template< typename LAI > @@ -301,36 +298,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 > @@ -358,7 +364,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 e4c1bdaa8a0..c093b556ad5 100644 --- a/src/coreComponents/linearAlgebra/interfaces/direct/SuperLUDist.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/direct/SuperLUDist.cpp @@ -22,6 +22,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" @@ -113,7 +114,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 @@ -130,7 +130,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 ); @@ -212,6 +211,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() ); @@ -221,10 +222,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; @@ -234,8 +235,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, @@ -248,8 +249,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 > @@ -320,7 +320,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 ); @@ -328,7 +328,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 ); @@ -424,7 +424,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 9f8f33bc659..bb00146dab6 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreExport.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreExport.cpp @@ -158,10 +158,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] ); } ); } @@ -176,7 +175,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/HypreInterface.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreInterface.cpp index 783c2e6b483..f9e5c42b89a 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreInterface.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreInterface.cpp @@ -109,7 +109,7 @@ geos::HypreInterface::createPreconditioner( LinearSolverParameters params ) std::unique_ptr< PreconditionerBase< HypreInterface > > geos::HypreInterface::createPreconditioner( LinearSolverParameters params, - array1d< HypreVector > const & nearNullKernel ) + arrayView1d< HypreVector const > nearNullKernel ) { return std::make_unique< HyprePreconditioner >( std::move( params ), nearNullKernel ); } diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreInterface.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreInterface.hpp index d6bdfd44d0a..382309b9f05 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreInterface.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreInterface.hpp @@ -70,7 +70,7 @@ struct HypreInterface */ static std::unique_ptr< PreconditionerBase< HypreInterface > > createPreconditioner( LinearSolverParameters params, - array1d< HypreVector > const & nearNullKernel ); + arrayView1d< HypreVector const > nearNullKernel ); /// Alias for HypreMatrix using ParallelMatrix = HypreMatrix; diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreKernels.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreKernels.hpp index b71248aa0d5..e4ba782d79e 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreKernels.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreKernels.hpp @@ -210,72 +210,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 a5342efdc7b..04e6408a4bd 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.cpp @@ -22,6 +22,7 @@ #include "common/TimingMacros.hpp" #include "common/GeosxConfig.hpp" #include "codingUtilities/Utilities.hpp" +#include "linearAlgebra/DofManager.hpp" #include "linearAlgebra/interfaces/hypre/HypreKernels.hpp" #include "linearAlgebra/interfaces/hypre/HypreUtils.hpp" #include "LvArray/src/output.hpp" @@ -690,17 +691,19 @@ void HypreMatrix::multiplyRAP( HypreMatrix const & R, dst.parCSRtoIJ( dst_parcsr ); } -void HypreMatrix::multiplyPtAP( HypreMatrix const & P, - HypreMatrix & dst ) const +void HypreMatrix::multiplyP1tAP2( 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 ); @@ -862,15 +865,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(); @@ -929,24 +929,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: { @@ -998,7 +988,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() ); @@ -1011,7 +1001,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 @@ -1055,6 +1045,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 @@ -1091,6 +1091,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 { @@ -1339,17 +1422,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 @@ -1365,13 +1452,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 7a268dfda03..12e7df67d3d 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMatrix.hpp @@ -21,7 +21,6 @@ #define GEOS_LINEARALGEBRA_INTERFACES_HYPREMATRIX_HPP_ #include "common/DataTypes.hpp" -#include "linearAlgebra/DofManager.hpp" #include "linearAlgebra/interfaces/hypre/HypreVector.hpp" #include "linearAlgebra/interfaces/hypre/HypreExport.hpp" #include "linearAlgebra/common/LinearOperator.hpp" @@ -130,6 +129,10 @@ 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; + using MatrixBase::multiplyP1tAP2; virtual void create( CRSMatrixView< real64 const, globalIndex const > const & localMatrix, localIndex const numLocalColumns, @@ -252,8 +255,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 multiplyP1tAP2( Matrix const & P1, + Matrix const & P2, + Matrix & dst ) const override; virtual void gemv( real64 const alpha, HypreVector const & x, @@ -303,20 +307,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 06d2a98981c..cc16aaaf59b 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.cpp @@ -45,7 +45,7 @@ struct HypreNullSpace namespace { -void convertRigidBodyModes( arrayView1d< HypreVector > const & nearNullKernel, +void convertRigidBodyModes( arrayView1d< HypreVector const > nearNullKernel, array1d< HYPRE_ParVector > & nullSpacePointer ) { if( nearNullKernel.empty() ) @@ -86,11 +86,14 @@ 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 ) ); // 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::iluk || 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 { @@ -250,6 +255,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; @@ -264,6 +275,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 ) @@ -273,7 +293,7 @@ HyprePreconditioner::HyprePreconditioner( LinearSolverParameters params ) {} HyprePreconditioner::HyprePreconditioner( LinearSolverParameters params, - arrayView1d< HypreVector > const & nearNullKernel ) + arrayView1d< HypreVector const > nearNullKernel ) : HyprePreconditioner( std::move( params ) ) { if( m_params.preconditionerType == LinearSolverParameters::PreconditionerType::amg && @@ -303,12 +323,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 ); diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.hpp index f26f15a636f..bb16b994797 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HyprePreconditioner.hpp @@ -63,7 +63,7 @@ class HyprePreconditioner final : public PreconditionerBase< HypreInterface > * @param nearNullKernel the user-provided near null kernel */ HyprePreconditioner( LinearSolverParameters params, - arrayView1d< HypreVector > const & nearNullKernel ); + arrayView1d< HypreVector const > nearNullKernel ); /** * @brief Destructor. diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreUtils.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreUtils.cpp index 80a6ffba199..2fbf79c7aff 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreUtils.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreUtils.cpp @@ -24,6 +24,8 @@ #include <_hypre_parcsr_mv.h> #include <_hypre_parcsr_ls.h> +#include + namespace geos { @@ -46,6 +48,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, @@ -146,7 +210,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; } @@ -165,6 +229,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 ); @@ -173,6 +238,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 47f901330f0..8db78a16c4e 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreUtils.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreUtils.hpp @@ -202,15 +202,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) @@ -284,6 +294,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 @@ -375,7 +429,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::iluk, 0 }, { LinearSolverParameters::AMG::SmootherType::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 048a299c3b7..56b3009022c 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/PetscExport.hpp b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscExport.hpp index e32eab2ab40..fc3a885f3a4 100644 --- a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscExport.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscExport.hpp @@ -26,7 +26,7 @@ extern "C" struct _p_IS; /// VecScatter struct forward declaration -extern "C" struct _p_VecScatter; +extern "C" struct _p_PetscSF; namespace geos { @@ -119,7 +119,7 @@ class PetscExport using IS = struct _p_IS *; /// Alias for PETSc vector scatter struct pointer - using VecScatter = struct _p_VecScatter *; + using VecScatter = struct _p_PetscSF *; /// Target rank for single-rank export/import integer m_targetRank = -1; diff --git a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscInterface.cpp b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscInterface.cpp index bb404c744eb..30671b4a040 100644 --- a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscInterface.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscInterface.cpp @@ -70,7 +70,7 @@ PetscInterface::createPreconditioner( LinearSolverParameters params ) std::unique_ptr< PreconditionerBase< PetscInterface > > PetscInterface::createPreconditioner( LinearSolverParameters params, - array1d< PetscVector > const & nearNullKernel ) + arrayView1d< PetscVector const > nearNullKernel ) { return std::make_unique< PetscPreconditioner >( params, nearNullKernel ); } diff --git a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscInterface.hpp b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscInterface.hpp index d49caa38cee..2767b7d5efb 100644 --- a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscInterface.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscInterface.hpp @@ -70,7 +70,7 @@ struct PetscInterface */ static std::unique_ptr< PreconditionerBase< PetscInterface > > createPreconditioner( LinearSolverParameters params, - array1d< PetscVector > const & nearNullKernel ); + arrayView1d< PetscVector const > nearNullKernel ); /// Alias for PetscMatrix using ParallelMatrix = PetscMatrix; diff --git a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.cpp b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.cpp index dac9a208a7c..32fa3d1a7f7 100644 --- a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.cpp @@ -24,6 +24,7 @@ #endif #include "codingUtilities/Utilities.hpp" +#include "linearAlgebra/DofManager.hpp" #include "linearAlgebra/interfaces/petsc/PetscUtils.hpp" #include "common/MpiWrapper.hpp" @@ -36,6 +37,8 @@ namespace geos { +static_assert( sizeof(globalIndex) == sizeof(PetscInt), "Mismatch in global index type width between GEOS and Petsc" ); + PetscMatrix::PetscMatrix() : LinearOperator(), MatrixBase() @@ -304,7 +307,7 @@ void PetscMatrix::close() // Clear any rows that have been marked for clearing if( assembled() ) { - GEOS_LAI_CHECK_ERROR( MatZeroRows( m_mat, m_rowsToClear.size(), m_rowsToClear.data(), 0.0, nullptr, nullptr ) ); + GEOS_LAI_CHECK_ERROR( MatZeroRows( m_mat, m_rowsToClear.size(), reinterpret_cast< PetscInt const * >(m_rowsToClear.data()), 0.0, nullptr, nullptr ) ); for( localIndex i = 0; i < m_rowsToClear.size(); ++i ) { @@ -655,8 +658,70 @@ void PetscMatrix::leftRightScale( PetscVector const & vecLeft, void PetscMatrix::computeScalingVector( PetscVector & scaling ) const { - GEOS_UNUSED_VAR( scaling ); - GEOS_ERROR( "Not implemented!!!" ); + GEOS_LAI_ASSERT( ready() ); + + // Get number of components + integer const numComp = m_dofManager->numComponents(); + + // Get local dof component labels + array1d< integer > dofLabels( numLocalRows() ); + m_dofManager->getLocalDofComponentLabels( dofLabels ); + + array1d< real64 > weights( numComp ); + + PetscInt firstRow, lastRow; + GEOS_LAI_CHECK_ERROR( MatGetOwnershipRange( m_mat, &firstRow, &lastRow ) ); + + PetscInt firstCol, lastCol; + GEOS_LAI_CHECK_ERROR( MatGetOwnershipRangeColumn( m_mat, &firstCol, &lastCol ) ); + + 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 ); + integer const rowLabel = dofLabels[localRow]; + + for( int k = 0; k < numEntries; ++k ) + { + PetscInt const globalCol = cols[k]; + if( firstCol <= globalCol && globalCol < lastCol ) + { + localIndex const localCol = LvArray::integerConversion< localIndex >( globalCol - firstCol ); + integer const colLabel = dofLabels[localCol]; + if( rowLabel == colLabel ) + { + weights[rowLabel] += vals[k]*vals[k]; + } + } + } + GEOS_LAI_CHECK_ERROR( MatRestoreRow( m_mat, row, &numEntries, &cols, &vals ) ); + } + + // Compute total squared Frobenius norm across all ranks + MpiWrapper::sum( Span< real64 const >{weights}, Span< real64 >{weights}, comm() ); + + // Compute the scaling weights + for( integer c = 0; c < numComp; ++c ) + { + // TODO: consider digit truncation like hypre does + weights[c] = LvArray::math::sqrt( LvArray::math::invSqrt( weights[c] ) ); + } + + // Populate the scaling vector + scaling.create( numLocalRows(), comm() ); + Vec & vec = scaling.unwrapped(); + PetscScalar *values; + GEOS_LAI_CHECK_ERROR( VecGetArray( vec, &values ) ); + for( PetscInt row = firstRow; row < lastRow; ++row ) + { + localIndex const localRow = LvArray::integerConversion< localIndex >( row - firstRow ); + values[localRow] = weights[dofLabels[localRow]]; + } + ; + GEOS_LAI_CHECK_ERROR( VecRestoreArray( vec, &values ) ); } void PetscMatrix::multiplyRAP( PetscMatrix const & R, @@ -666,13 +731,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; } @@ -681,18 +745,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; } @@ -708,17 +768,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 ); }; @@ -737,7 +796,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_GEOS ); @@ -792,7 +851,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; @@ -858,7 +917,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 ); @@ -867,16 +926,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; } @@ -891,6 +951,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 @@ -899,17 +982,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 @@ -922,6 +1007,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, reinterpret_cast< globalIndex const * >(&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 { @@ -931,11 +1108,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; } diff --git a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.hpp b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.hpp index 847ed89cc0a..7add0fa7bc6 100644 --- a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscMatrix.hpp @@ -119,6 +119,10 @@ 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; + using MatrixBase::multiplyP1tAP2; virtual void createWithLocalSize( localIndex const localRows, localIndex const localCols, @@ -283,14 +287,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 822bf336aae..09d9138b48c 100644 --- a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscPreconditioner.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscPreconditioner.cpp @@ -30,7 +30,7 @@ namespace { void convertRigidBodyModes( LinearSolverParameters const & params, - array1d< PetscVector > const & nearNullKernel, + arrayView1d< PetscVector const > nearNullKernel, MatNullSpace & nullsp ) { if( nearNullKernel.empty() ) @@ -62,7 +62,7 @@ PCType getPetscSmootherType( LinearSolverParameters::PreconditionerType const & static std::map< LinearSolverParameters::PreconditionerType, PCType > const typeMap = { { LinearSolverParameters::PreconditionerType::iluk, PCILU }, - { LinearSolverParameters::PreconditionerType::ic, PCICC }, + { LinearSolverParameters::PreconditionerType::ick, PCICC }, { LinearSolverParameters::PreconditionerType::jacobi, PCJACOBI }, { LinearSolverParameters::PreconditionerType::l1jacobi, PCJACOBI }, { LinearSolverParameters::PreconditionerType::fgs, PCSOR }, @@ -144,14 +144,6 @@ void createPetscAMG( LinearSolverParameters const & params, // Default options only for the moment GEOS_LAI_CHECK_ERROR( PCSetType( precond, PCGAMG ) ); - if( !params.isSymmetric ) - { - // Usually GEOSX matrix is not symmetric, but GAMG is designed for symmetric matrices - // In case of a general matrix, we need to compute a symmetric graph (slightly heavier - // than the default, but necessary in case of asymmetric matrices). - GEOS_LAI_CHECK_ERROR( PCGAMGSetSymGraph( precond, PETSC_TRUE ) ); - } - // Add user-defined null space / rigid body mode support if( params.amg.nullSpaceType == LinearSolverParameters::AMG::NullSpaceType::rigidBodyModes && nullsp ) { @@ -167,6 +159,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 ) ); @@ -251,7 +246,7 @@ PetscPreconditioner::PetscPreconditioner( LinearSolverParameters params ) { } PetscPreconditioner::PetscPreconditioner( LinearSolverParameters params, - array1d< Vector > const & nearNullKernel ) + arrayView1d< Vector const > nearNullKernel ) : Base{}, m_params( std::move( params ) ), m_precond{}, @@ -327,7 +322,7 @@ void PetscPreconditioner::setup( PetscMatrix const & mat ) case LinearSolverParameters::PreconditionerType::bgs: case LinearSolverParameters::PreconditionerType::sgs: case LinearSolverParameters::PreconditionerType::iluk: - case LinearSolverParameters::PreconditionerType::ic: + case LinearSolverParameters::PreconditionerType::ick: { createPetscSmoother( m_params, m_precond ); break; diff --git a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscPreconditioner.hpp b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscPreconditioner.hpp index ac4a9891247..0d93db926a4 100644 --- a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscPreconditioner.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscPreconditioner.hpp @@ -68,7 +68,8 @@ class PetscPreconditioner final : public PreconditionerBase< PetscInterface > * @param params preconditioner parameters * @param nearNullKernel the user-provided near null kernel */ - PetscPreconditioner( LinearSolverParameters params, array1d< Vector > const & nearNullKernel ); + PetscPreconditioner( LinearSolverParameters params, + arrayView1d< Vector const > nearNullKernel ); /** * @brief Destructor. diff --git a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscSolver.cpp b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscSolver.cpp index 1470cd2e4a9..a8fbcdff3e1 100644 --- a/src/coreComponents/linearAlgebra/interfaces/petsc/PetscSolver.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/petsc/PetscSolver.cpp @@ -95,7 +95,7 @@ void PetscSolver::setup( PetscMatrix const & mat ) { // cast needed because of "void *" vs "PetscViewerAndFormat *" in last parameter using MonitorFunc = PetscErrorCode ( * )( KSP, PetscInt, PetscReal, void * ); - GEOS_LAI_CHECK_ERROR( KSPMonitorSet( m_solver, ( MonitorFunc ) KSPMonitorDefault, nullptr, nullptr ) ); + GEOS_LAI_CHECK_ERROR( KSPMonitorSet( m_solver, ( MonitorFunc ) KSPMonitorResidual, nullptr, nullptr ) ); } // This can be used to extract residual norm history: diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.cpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.cpp index 295714afa91..77f2d27bc9a 100644 --- a/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.cpp @@ -21,6 +21,8 @@ #include "codingUtilities/RTTypes.hpp" #include "codingUtilities/Utilities.hpp" +#include "common/GEOS_RAJA_Interface.hpp" +#include "linearAlgebra/DofManager.hpp" #include "linearAlgebra/interfaces/trilinos/EpetraUtils.hpp" #include @@ -31,6 +33,7 @@ #include #include #include +#include #include @@ -455,8 +458,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 ) ); @@ -471,21 +474,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, @@ -533,8 +537,62 @@ void EpetraMatrix::leftRightScale( EpetraVector const & vecLeft, void EpetraMatrix::computeScalingVector( EpetraVector & scaling ) const { - GEOS_UNUSED_VAR( scaling ); - GEOS_ERROR( "Not implemented!!!" ); + GEOS_LAI_ASSERT( ready() ); + + // Get number of components + integer const numComp = m_dofManager->numComponents(); + + // Get local dof component labels + array1d< integer > dofLabels( numLocalRows() ); + m_dofManager->getLocalDofComponentLabels( dofLabels ); + + localIndex firstCol = jlower(); + localIndex lastCol = jupper(); + + array1d< real64 > weights( numComp ); + Epetra_CrsMatrix const & mat = unwrapped(); + + // Compute local diagonal block squared Frobeniums norm for each dof component + forAll< parallelHostPolicy >( mat.NumMyRows(), [=, &mat, &weights, &dofLabels]( int const localRow ) + { + integer const rowLabel = dofLabels[localRow]; + int numEntries; + int * columns; + double * values; + mat.ExtractMyRowView( localRow, numEntries, values, columns ); + for( int k = 0; k < numEntries; ++k ) + { + long long const globalCol = mat.GCID64( columns[k] ); + if( firstCol <= globalCol && globalCol < lastCol ) + { + localIndex const localCol = LvArray::integerConversion< localIndex >( globalCol - firstCol ); + integer const colLabel = dofLabels[localCol]; + if( rowLabel == colLabel ) + { + RAJA::atomicAdd( parallelHostAtomic{}, &weights[rowLabel], values[k] * values[k] ); + } + } + } + } ); + + // Compute total squared Frobenius norm across all ranks + MpiWrapper::sum( Span< real64 const >{weights}, Span< real64 >{weights}, comm() ); + + // Compute the scaling weights + for( integer c = 0; c < numComp; ++c ) + { + // TODO: consider digit truncation like hypre does + weights[c] = LvArray::math::sqrt( LvArray::math::invSqrt( weights[c] ) ); + } + + // Populate the scaling vector + scaling.create( numLocalRows(), comm() ); + Epetra_Vector & vec = scaling.unwrapped(); + real64 * const values = vec.Values(); + forAll< parallelHostPolicy >( vec.MyLength(), [=, &mat, &weights, &dofLabels]( int const localRow ) + { + values[localRow] = weights[dofLabels[localRow]]; + } ); } void EpetraMatrix::transpose( EpetraMatrix & dst ) const @@ -551,21 +609,19 @@ void EpetraMatrix::transpose( EpetraMatrix & dst ) const void EpetraMatrix::separateComponentFilter( EpetraMatrix & dst, integer const dofsPerNode ) const { - localIndex const maxRowEntries = maxRowLength(); - integer const remainder = maxRowEntries % dofsPerNode; - GEOS_LAI_ASSERT_EQ( remainder, 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; @@ -583,7 +639,7 @@ void EpetraMatrix::separateComponentFilter( EpetraMatrix & dst, } } ); - dst.create( tempMatView.toViewConst(), numLocalCols(), MPI_COMM_GEOS ); + dst.create( tempMat.toViewConst(), numLocalCols(), MPI_COMM_GEOS ); dst.setDofManager( dofManager() ); } @@ -624,39 +680,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 @@ -672,7 +780,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 ) ); @@ -734,6 +842,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() ); @@ -749,9 +863,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] = 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 ) { - lengths[localRow] = m_matrix->NumMyEntries( LvArray::integerConversion< int >( localRow ) ); + 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 ); } ); } @@ -773,7 +908,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() ); } @@ -787,6 +922,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 { @@ -824,6 +1032,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 49b30321dfe..4abecdab9d3 100644 --- a/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/EpetraMatrix.hpp @@ -108,6 +108,10 @@ 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; + using MatrixBase::multiplyP1tAP2; virtual void createWithLocalSize( localIndex const localRows, localIndex const localCols, @@ -269,6 +273,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 */ @@ -278,12 +287,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/TrilinosInterface.cpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosInterface.cpp index 39c5883876f..cb40407aeb4 100644 --- a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosInterface.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosInterface.cpp @@ -61,7 +61,7 @@ TrilinosInterface::createPreconditioner( LinearSolverParameters params ) std::unique_ptr< PreconditionerBase< TrilinosInterface > > TrilinosInterface::createPreconditioner( LinearSolverParameters params, - array1d< EpetraVector > const & nearNullKernel ) + arrayView1d< EpetraVector const > nearNullKernel ) { return std::make_unique< TrilinosPreconditioner >( params, nearNullKernel ); } diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosInterface.hpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosInterface.hpp index 6c3a2a06925..d041af86836 100644 --- a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosInterface.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosInterface.hpp @@ -70,7 +70,7 @@ struct TrilinosInterface */ static std::unique_ptr< PreconditionerBase< TrilinosInterface > > createPreconditioner( LinearSolverParameters params, - array1d< EpetraVector > const & nearNullKernel ); + arrayView1d< EpetraVector const > nearNullKernel ); /// Alias for EpetraMatrix using ParallelMatrix = EpetraMatrix; diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosPreconditioner.cpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosPreconditioner.cpp index 1f2bfb02faf..6d5f48dccc6 100644 --- a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosPreconditioner.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosPreconditioner.cpp @@ -37,7 +37,7 @@ TrilinosPreconditioner::TrilinosPreconditioner( LinearSolverParameters params ) m_nullSpacePointer{} {} -void convertRigidBodyModes( array1d< EpetraVector > const & nearNullKernel, +void convertRigidBodyModes( arrayView1d< EpetraVector const > nearNullKernel, array2d< real64 > & nullSpacePointer ) { if( nearNullKernel.empty() ) @@ -62,7 +62,7 @@ void convertRigidBodyModes( array1d< EpetraVector > const & nearNullKernel, } TrilinosPreconditioner::TrilinosPreconditioner( LinearSolverParameters params, - array1d< Vector > const & nearNullKernel ) + arrayView1d< Vector const > nearNullKernel ) : Base{}, m_params( std::move( params ) ), m_precond{}, @@ -104,8 +104,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::ick, "IC" }, + { LinearSolverParameters::AMG::SmootherType::iluk, "ILU" }, { LinearSolverParameters::AMG::SmootherType::ilut, "ILUT" }, }; @@ -155,14 +155,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 ) { @@ -176,10 +179,7 @@ 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 ) @@ -188,13 +188,15 @@ Ifpack::EPrecType getIfpackPrecondType( LinearSolverParameters::PreconditionerTy { { LinearSolverParameters::PreconditionerType::iluk, Ifpack::ILU }, { LinearSolverParameters::PreconditionerType::ilut, Ifpack::ILUT }, - { LinearSolverParameters::PreconditionerType::ic, Ifpack::IC }, + { LinearSolverParameters::PreconditionerType::ick, 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 } }; @@ -203,11 +205,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 ); @@ -246,6 +253,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; @@ -267,11 +279,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 ); @@ -287,13 +307,15 @@ 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::ilut: - case LinearSolverParameters::PreconditionerType::ic: + case LinearSolverParameters::PreconditionerType::ick: case LinearSolverParameters::PreconditionerType::ict: case LinearSolverParameters::PreconditionerType::direct: { @@ -302,7 +324,6 @@ void TrilinosPreconditioner::setup( Matrix const & mat ) } case LinearSolverParameters::PreconditionerType::none: { - m_precond.reset(); break; } default: diff --git a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosPreconditioner.hpp b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosPreconditioner.hpp index 2efdf01b192..ddaa81da614 100644 --- a/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosPreconditioner.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/trilinos/TrilinosPreconditioner.hpp @@ -59,7 +59,7 @@ class TrilinosPreconditioner final : public PreconditionerBase< TrilinosInterfac * @param nearNullKernel the user-provided near null kernel */ TrilinosPreconditioner( LinearSolverParameters params, - array1d< Vector > const & nearNullKernel ); + arrayView1d< Vector const > nearNullKernel ); /** * @brief Destructor. diff --git a/src/coreComponents/linearAlgebra/solvers/BicgstabSolver.cpp b/src/coreComponents/linearAlgebra/solvers/BicgstabSolver.cpp index e8e486916d8..13555a07b75 100644 --- a/src/coreComponents/linearAlgebra/solvers/BicgstabSolver.cpp +++ b/src/coreComponents/linearAlgebra/solvers/BicgstabSolver.cpp @@ -20,6 +20,7 @@ #include "BicgstabSolver.hpp" +#include "common/TimingMacros.hpp" #include "common/Stopwatch.hpp" #include "linearAlgebra/interfaces/InterfaceTypes.hpp" #include "common/LinearOperator.hpp" @@ -40,6 +41,7 @@ template< typename VECTOR > void BicgstabSolver< VECTOR >::solve( Vector const & b, Vector & x ) const { + GEOS_MARK_FUNCTION; Stopwatch watch; // Define vectors diff --git a/src/coreComponents/linearAlgebra/solvers/BicgstabSolver.hpp b/src/coreComponents/linearAlgebra/solvers/BicgstabSolver.hpp index 41400be9b27..0dd3d0a5390 100644 --- a/src/coreComponents/linearAlgebra/solvers/BicgstabSolver.hpp +++ b/src/coreComponents/linearAlgebra/solvers/BicgstabSolver.hpp @@ -35,7 +35,7 @@ namespace geos * from Y. Saad (2003). */ template< typename VECTOR > -class BicgstabSolver : public KrylovSolver< VECTOR > +class BicgstabSolver final : public KrylovSolver< VECTOR > { public: @@ -72,9 +72,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 d77d4a7dbfc..4d500926c0b 100644 --- a/src/coreComponents/linearAlgebra/solvers/BlockPreconditioner.cpp +++ b/src/coreComponents/linearAlgebra/solvers/BlockPreconditioner.cpp @@ -28,13 +28,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 }, @@ -43,24 +39,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 ); } } @@ -76,16 +74,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 ); @@ -105,14 +131,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(); @@ -125,7 +151,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 ) ); @@ -134,7 +160,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; @@ -152,9 +178,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 ); @@ -171,18 +194,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.multiplyP1tAP2( *m_prolongators[0], *m_prolongators[1], m_matBlocks( 0, 1 ) ); + mat.multiplyP1tAP2( *m_prolongators[1], *m_prolongators[0], m_matBlocks( 1, 0 ) ); } applyBlockScaling(); @@ -195,36 +225,44 @@ 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; - for( localIndex i = 0; i < 2; ++i ) + m_prolongators[0]->applyTranspose( src, m_rhs( 0 ) ); + m_prolongators[1]->applyTranspose( src, m_rhs( 1 ) ); + + if( m_params.scaling != LinearSolverParameters::Block::Scaling::None ) { - m_rhs( i ).scale( m_scaling[i] ); + 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 > @@ -233,8 +271,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 57b6fe7929b..58622d2f933 100644 --- a/src/coreComponents/linearAlgebra/solvers/BlockPreconditioner.hpp +++ b/src/coreComponents/linearAlgebra/solvers/BlockPreconditioner.hpp @@ -24,43 +24,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 @@ -113,18 +81,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, + stdVector< DofManager::SubComponent > blockDofs, + std::unique_ptr< PreconditionerBase< LAI > > solver, + real64 const scaling = 1.0 ); /** * @brief Setup data for one of the two blocks. @@ -133,22 +108,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, - stdVector< DofManager::SubComponent > blockDofs, - std::unique_ptr< PreconditionerBase< LAI > > solver, + std::vector< DofManager::SubComponent > blockDofs, + 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 @@ -168,14 +151,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). @@ -187,32 +177,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< stdVector< DofManager::SubComponent >, 2 > m_blockDofs; + std::array< stdVector< 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 ca80ddb1c5a..8a62ac39304 100644 --- a/src/coreComponents/linearAlgebra/solvers/CgSolver.cpp +++ b/src/coreComponents/linearAlgebra/solvers/CgSolver.cpp @@ -21,6 +21,7 @@ #include "CgSolver.hpp" +#include "common/TimingMacros.hpp" #include "common/Stopwatch.hpp" #include "linearAlgebra/interfaces/InterfaceTypes.hpp" #include "common/LinearOperator.hpp" @@ -58,6 +59,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 diff --git a/src/coreComponents/linearAlgebra/solvers/CgSolver.hpp b/src/coreComponents/linearAlgebra/solvers/CgSolver.hpp index bc88229735d..375c0e9f6af 100644 --- a/src/coreComponents/linearAlgebra/solvers/CgSolver.hpp +++ b/src/coreComponents/linearAlgebra/solvers/CgSolver.hpp @@ -35,7 +35,7 @@ namespace geos * from Y. Saad (2003). */ template< typename VECTOR > -class CgSolver : public KrylovSolver< VECTOR > +class CgSolver final : public KrylovSolver< VECTOR > { public: @@ -72,9 +72,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"; }; @@ -97,6 +97,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 f2b801c5d04..f3ce099b011 100644 --- a/src/coreComponents/linearAlgebra/solvers/GmresSolver.cpp +++ b/src/coreComponents/linearAlgebra/solvers/GmresSolver.cpp @@ -19,6 +19,7 @@ #include "GmresSolver.hpp" +#include "common/TimingMacros.hpp" #include "common/Stopwatch.hpp" #include "linearAlgebra/interfaces/InterfaceTypes.hpp" #include "linearAlgebra/solvers/KrylovUtils.hpp" @@ -31,11 +32,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 @@ -90,16 +90,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; @@ -127,6 +132,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 ) { @@ -167,6 +176,8 @@ void GmresSolver< VECTOR >::solve( Vector const & b, H( j+1, j ) = w.norm2(); 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 1f41810c7d5..855510a5b84 100644 --- a/src/coreComponents/linearAlgebra/solvers/GmresSolver.hpp +++ b/src/coreComponents/linearAlgebra/solvers/GmresSolver.hpp @@ -35,7 +35,7 @@ namespace geos * from Y. Saad (2003). */ template< typename VECTOR > -class GmresSolver : public KrylovSolver< VECTOR > +class GmresSolver final : public KrylovSolver< VECTOR > { public: @@ -72,9 +72,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"; }; @@ -96,10 +96,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 c2c7eb4d59d..a4c2f6c6570 100644 --- a/src/coreComponents/linearAlgebra/solvers/KrylovSolver.cpp +++ b/src/coreComponents/linearAlgebra/solvers/KrylovSolver.cpp @@ -21,6 +21,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 @@ -68,6 +69,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 ); @@ -82,8 +89,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 44856621892..dc50100a21c 100644 --- a/src/coreComponents/linearAlgebra/solvers/KrylovSolver.hpp +++ b/src/coreComponents/linearAlgebra/solvers/KrylovSolver.hpp @@ -13,6 +13,10 @@ * ------------------------------------------------------------------------------------------------------------ */ +/** + * @file KrylovSolver.hpp + */ + #ifndef GEOS_LINEARALGEBRA_SOLVERS_KRYLOVSOLVER_HPP_ #define GEOS_LINEARALGEBRA_SOLVERS_KRYLOVSOLVER_HPP_ @@ -47,9 +51,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. @@ -61,11 +66,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. @@ -73,10 +73,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/PreconditionerIdentity.hpp b/src/coreComponents/linearAlgebra/solvers/PreconditionerIdentity.hpp index 696a8ab481d..0fa06b4fc9b 100644 --- a/src/coreComponents/linearAlgebra/solvers/PreconditionerIdentity.hpp +++ b/src/coreComponents/linearAlgebra/solvers/PreconditionerIdentity.hpp @@ -13,10 +13,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 @@ -40,8 +43,6 @@ class PreconditionerIdentity : public PreconditionerBase< LAI > /// Alias for matrix type using Matrix = typename Base::Matrix; - virtual ~PreconditionerIdentity() = default; - /** * @brief Apply operator to a vector. * @@ -57,6 +58,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..2dce99c773d --- /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 GEOS_LINEARALGEBRA_SOLVERS_PRECONDITIONERNULL_HPP_ +#define GEOS_LINEARALGEBRA_SOLVERS_PRECONDITIONERNULL_HPP_ + +#include "linearAlgebra/common/PreconditionerBase.hpp" + +namespace geos +{ + +/** + * @brief Common interface for null (zero) 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 //GEOS_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..fb3c145ad57 --- /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 GEOS_USE_TRILINOS +template class RichardsonSolver< TrilinosInterface::ParallelVector >; +template class RichardsonSolver< BlockVectorView< TrilinosInterface::ParallelVector > >; +#endif + +#ifdef GEOS_USE_HYPRE +template class RichardsonSolver< HypreInterface::ParallelVector >; +template class RichardsonSolver< BlockVectorView< HypreInterface::ParallelVector > >; +#endif + +#ifdef GEOS_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..14db37190fc --- /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 GEOS_LINEARALGEBRA_SOLVERS_RICHARDSONSOLVER_HPP_ +#define GEOS_LINEARALGEBRA_SOLVERS_RICHARDSONSOLVER_HPP_ + +#include "linearAlgebra/solvers/KrylovSolver.hpp" + +namespace geos +{ + +/** + * @brief Implements 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 //GEOS_LINEARALGEBRA_SOLVERS_RICHARDSONSOLVER_HPP_ diff --git a/src/coreComponents/linearAlgebra/unitTests/testKrylovSolvers.cpp b/src/coreComponents/linearAlgebra/unitTests/testKrylovSolvers.cpp index 898cfa63f41..90b4a4cdab0 100644 --- a/src/coreComponents/linearAlgebra/unitTests/testKrylovSolvers.cpp +++ b/src/coreComponents/linearAlgebra/unitTests/testKrylovSolvers.cpp @@ -56,6 +56,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 { @@ -155,10 +165,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 GEOS_USE_TRILINOS INSTANTIATE_TYPED_TEST_SUITE_P( Trilinos, KrylovSolverTest, TrilinosInterface, ); @@ -245,10 +261,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 GEOS_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 0154435afe1..d9404b5e7b4 100644 --- a/src/coreComponents/linearAlgebra/unitTests/testLinearSolverParametersEnums.cpp +++ b/src/coreComponents/linearAlgebra/unitTests/testLinearSolverParametersEnums.cpp @@ -34,6 +34,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 ) ); } @@ -51,7 +52,7 @@ TEST( LinearSolverParametersEnums, PreconditionerType ) ASSERT_EQ( "chebyshev", toString( EnumType::chebyshev ) ); ASSERT_EQ( "iluk", toString( EnumType::iluk ) ); ASSERT_EQ( "ilut", toString( EnumType::ilut ) ); - ASSERT_EQ( "icc", toString( EnumType::ic ) ); // Notice the discrepancy here + ASSERT_EQ( "ick", toString( EnumType::ick ) ); // Notice the discrepancy here ASSERT_EQ( "ict", toString( EnumType::ict ) ); ASSERT_EQ( "amg", toString( EnumType::amg ) ); ASSERT_EQ( "mgr", toString( EnumType::mgr ) ); @@ -135,9 +136,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( "iluk", toString( EnumType::iluk ) ); ASSERT_EQ( "ilut", toString( EnumType::ilut ) ); - ASSERT_EQ( "ic0", toString( EnumType::ic0 ) ); + ASSERT_EQ( "ick", toString( EnumType::ick ) ); ASSERT_EQ( "ict", toString( EnumType::ict ) ); } diff --git a/src/coreComponents/linearAlgebra/utilities/BlockOperator.hpp b/src/coreComponents/linearAlgebra/utilities/BlockOperator.hpp index f25068efeb4..75dc85db741 100644 --- a/src/coreComponents/linearAlgebra/utilities/BlockOperator.hpp +++ b/src/coreComponents/linearAlgebra/utilities/BlockOperator.hpp @@ -58,17 +58,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(); @@ -107,14 +96,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 a7da7ca9a79..f3620e0048a 100644 --- a/src/coreComponents/linearAlgebra/utilities/BlockOperatorView.hpp +++ b/src/coreComponents/linearAlgebra/utilities/BlockOperatorView.hpp @@ -55,11 +55,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 977135fb419..48f5663931b 100644 --- a/src/coreComponents/linearAlgebra/utilities/InverseNormalOperator.hpp +++ b/src/coreComponents/linearAlgebra/utilities/InverseNormalOperator.hpp @@ -53,7 +53,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 fdd0c8e323a..0e014d430b8 100644 --- a/src/coreComponents/linearAlgebra/utilities/LAIHelperFunctions.hpp +++ b/src/coreComponents/linearAlgebra/utilities/LAIHelperFunctions.hpp @@ -206,144 +206,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, - stdVector< 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_GEOS ); + rigidBodyModes[k].create( numLocalDof, MPI_COMM_GEOS ); 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_GEOS ); - { - 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_GEOS ); + 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_GEOS ); - { - 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_GEOS ); - { - 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_GEOS ); - { - 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 07d95955bea..1507462f31c 100644 --- a/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp +++ b/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp @@ -21,6 +21,7 @@ #define GEOS_LINEARALGEBRA_UTILITIES_LINEARSOLVERPARAMETERS_HPP_ #include "common/format/EnumStrings.hpp" +#include "common/DataTypes.hpp" namespace geos { @@ -43,6 +44,7 @@ struct LinearSolverParameters gmres, ///< GMRES fgmres, ///< Flexible GMRES bicgstab, ///< BiCGStab + richardson, ///< Richardson iteration preconditioner ///< Preconditioner only }; @@ -60,13 +62,13 @@ struct LinearSolverParameters chebyshev, ///< Chebyshev polynomial smoothing iluk, ///< Incomplete LU with k-level of fill ilut, ///< Incomplete LU with thresholding - ic, ///< Incomplete Cholesky + ick, ///< Incomplete Cholesky with k-level of fill ict, ///< Incomplete Cholesky with thresholding amg, ///< Algebraic Multigrid 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] @@ -130,6 +132,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 { @@ -167,9 +184,9 @@ struct LinearSolverParameters sgs, ///< Symmetric Gauss-Seidel smoothing l1sgs, ///< l1-Symmetric Gauss-Seidel smoothing chebyshev, ///< Chebyshev polynomial smoothing - ilu0, ///< ILU(0) + iluk, ///< Incomplete LU with k-level of fill ilut, ///< Incomplete LU with thresholding - ic0, ///< Incomplete Cholesky + ick, ///< Incomplete Cholesky with k-level of fill ict ///< Incomplete Cholesky with thresholding }; @@ -246,6 +263,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 @@ -255,6 +273,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. integer aggressiveInterpMaxNonZeros = 16; ///< Aggressive Interpolation - Max. nonzeros/row. PreOrPost preOrPostSmoothing = PreOrPost::both; ///< Pre and/or post smoothing @@ -287,7 +306,7 @@ struct LinearSolverParameters compositionalMultiphaseHybridFVM, ///< hybrid finite volume compositional multiphase flow compositionalMultiphaseReservoirFVM, ///< finite volume compositional multiphase flow with wells compositionalMultiphaseReservoirHybridFVM, ///< hybrid finite volume compositional multiphase flow with wells - immiscibleMultiphaseFVM, ///< finite volume immiscible multiphase flow + immiscibleMultiphaseFVM, ///< finite volume immiscible multiphase flow reactiveCompositionalMultiphaseOBL, ///< finite volume reactive compositional flow with OBL thermalCompositionalMultiphaseFVM, ///< finite volume thermal compositional multiphase flow thermalCompositionalMultiphaseReservoirFVM,///< finite volume thermal compositional multiphase flow @@ -303,8 +322,8 @@ 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 - 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 that jacobi can be applied to the well + ///< block } mgr; ///< Multigrid reduction (MGR) parameters @@ -314,7 +333,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 @@ -322,6 +341,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. @@ -331,6 +402,7 @@ ENUM_STRINGS( LinearSolverParameters::SolverType, "gmres", "fgmres", "bicgstab", + "richardson", "preconditioner" ); /// Declare strings associated with enumeration values. @@ -344,7 +416,7 @@ ENUM_STRINGS( LinearSolverParameters::PreconditionerType, "chebyshev", "iluk", "ilut", - "icc", + "ick", "ict", "amg", "mgr", @@ -416,9 +488,9 @@ ENUM_STRINGS( LinearSolverParameters::AMG::SmootherType, "sgs", "l1sgs", "chebyshev", - "ilu0", + "iluk", "ilut", - "ic0", + "ick", "ict" ); /// Declare strings associated with enumeration values. @@ -475,6 +547,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 2382af2e0c2..96a53529eb8 100644 --- a/src/coreComponents/linearAlgebra/utilities/NormalOperator.hpp +++ b/src/coreComponents/linearAlgebra/utilities/NormalOperator.hpp @@ -16,8 +16,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" @@ -28,19 +28,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 @@ -50,11 +50,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 @@ -116,4 +111,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 18d1ed2c35a..cd30c6e08f0 100644 --- a/src/coreComponents/linearAlgebra/utilities/TransposeOperator.hpp +++ b/src/coreComponents/linearAlgebra/utilities/TransposeOperator.hpp @@ -17,8 +17,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" @@ -27,21 +27,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. @@ -52,11 +52,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) @@ -70,8 +65,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 { @@ -79,8 +73,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 { @@ -88,8 +81,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 { @@ -97,8 +89,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 { @@ -121,4 +112,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 433d55155ef..372eb398618 100644 --- a/src/coreComponents/mesh/CMakeLists.txt +++ b/src/coreComponents/mesh/CMakeLists.txt @@ -86,7 +86,6 @@ set( mesh_headers generators/InternalWellboreGenerator.hpp generators/MeshComponentBase.hpp generators/MeshGeneratorBase.hpp - generators/ParMETISInterface.hpp generators/ParticleMeshGenerator.hpp generators/PartitionDescriptor.hpp generators/PrismUtilities.hpp @@ -165,7 +164,6 @@ set( mesh_sources generators/InternalWellboreGenerator.cpp generators/MeshComponentBase.cpp generators/MeshGeneratorBase.cpp - generators/ParMETISInterface.cpp generators/ParticleMeshGenerator.cpp generators/Region.cpp generators/WellGeneratorBase.cpp @@ -186,7 +184,8 @@ set( mesh_sources simpleGeometricObjects/ThickPlane.cpp utilities/ComputationalGeometry.cpp ) -set( dependencyList ${parallelDeps} schema constitutive finiteElement parmetis metis ) +set( dependencyList ${parallelDeps} schema constitutive finiteElement ) +set( tplDependencyList "" ) if( ENABLE_VTK ) message(STATUS "Adding VTK readers") @@ -208,16 +207,22 @@ if( ENABLE_VTK ) generators/VTKWellGenerator.cpp generators/VTKUtilities.cpp ) - list( APPEND dependencyList VTK::IOLegacy VTK::FiltersParallelDIY2 ) + list( APPEND tplDependencyList VTK::IOLegacy VTK::FiltersParallelDIY2 ) if( ENABLE_MPI ) - list( APPEND dependencyList VTK::IOParallelXML VTK::ParallelMPI ) + list( APPEND tplDependencyList VTK::IOParallelXML VTK::ParallelMPI ) endif() endif() +if( ENABLE_PARMETIS ) + set( mesh_headers ${mesh_headers} generators/ParMETISInterface.hpp ) + set( mesh_sources ${mesh_sources} generators/ParMETISInterface.cpp ) + list( APPEND tplDependencyList parmetis metis ) +endif() + if( ENABLE_SCOTCH ) set( mesh_headers ${mesh_headers} generators/PTScotchInterface.hpp ) set( mesh_sources ${mesh_sources} generators/PTScotchInterface.cpp ) - list( APPEND dependencyList ptscotch ) + list( APPEND tplDependencyList SCOTCH::ptscotch SCOTCH::scotcherr ) endif() geos_decorate_link_dependencies( LIST decoratedDependencies @@ -226,7 +231,7 @@ geos_decorate_link_dependencies( LIST decoratedDependencies blt_add_library( NAME mesh SOURCES ${mesh_sources} HEADERS ${mesh_headers} - DEPENDS_ON ${decoratedDependencies} + DEPENDS_ON ${decoratedDependencies} ${tplDependencyList} OBJECT ${GEOS_BUILD_OBJ_LIBS} SHARED ${GEOS_BUILD_SHARED_LIBS} ) diff --git a/src/coreComponents/mesh/DomainPartition.hpp b/src/coreComponents/mesh/DomainPartition.hpp index dd108c46cd8..6644915738c 100644 --- a/src/coreComponents/mesh/DomainPartition.hpp +++ b/src/coreComponents/mesh/DomainPartition.hpp @@ -274,6 +274,21 @@ class DomainPartition : public dataRepository::Group stdVector< 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 349a4c91a31..f8692a95664 100644 --- a/src/coreComponents/mesh/ElementRegionManager.cpp +++ b/src/coreComponents/mesh/ElementRegionManager.cpp @@ -64,8 +64,7 @@ void ElementRegionManager::setMaxGlobalIndex() { m_localMaxGlobalIndex = std::max( m_localMaxGlobalIndex, subRegion.maxGlobalIndex() ); } ); - - m_maxGlobalIndex = MpiWrapper::max( m_localMaxGlobalIndex, MPI_COMM_GEOS ); + ObjectManagerBase::setMaxGlobalIndex(); } auto const & getUserAvailableKeys() diff --git a/src/coreComponents/physicsSolvers/LinearSolverParameters.cpp b/src/coreComponents/physicsSolvers/LinearSolverParameters.cpp index f5a17593011..09f07c089e1 100644 --- a/src/coreComponents/physicsSolvers/LinearSolverParameters.cpp +++ b/src/coreComponents/physicsSolvers/LinearSolverParameters.cpp @@ -23,8 +23,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 ) : @@ -36,13 +107,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 ). @@ -63,13 +134,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 ). @@ -130,6 +201,26 @@ LinearSolverParametersInput::LinearSolverParametersInput( string const & name, setInputFlag( InputFlags::OPTIONAL ). setDescription( "Exponent parameter 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 ). @@ -139,7 +230,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 ). @@ -150,7 +241,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 ). @@ -190,6 +281,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 ). @@ -203,8 +299,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 ). @@ -216,6 +312,8 @@ LinearSolverParametersInput::LinearSolverParametersInput( string const & name, setInputFlag( InputFlags::OPTIONAL ). setDescription( "ILU(T) threshold factor" ); + registerInputBlock< BlockParametersInput >( this, groupKeyStruct::blockString(), m_parameters.block ); + addLogLevel< logInfo::LinearSolver >(); addLogLevel< logInfo::LinearSolverConfiguration >(); } @@ -279,7 +377,16 @@ void LinearSolverParametersInput::postInputInitialization() // TODO input validation for other AMG parameters ? if( getLogLevel() > 0 ) + { print(); + } +} + +Group * LinearSolverParametersInput::createChild( string const & childKey, + string const & childName ) +{ + GEOS_UNUSED_VAR( childKey, childName ); + return nullptr; } void LinearSolverParametersInput::print() diff --git a/src/coreComponents/physicsSolvers/LinearSolverParameters.hpp b/src/coreComponents/physicsSolvers/LinearSolverParameters.hpp index 4b1699a1615..4d61ff97133 100644 --- a/src/coreComponents/physicsSolvers/LinearSolverParameters.hpp +++ b/src/coreComponents/physicsSolvers/LinearSolverParameters.hpp @@ -43,17 +43,9 @@ 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"; } @@ -62,6 +54,8 @@ class LinearSolverParametersInput : public dataRepository::Group void print(); + virtual Group * createChild( string const & childKey, string const & childName ) override final; + LinearSolverParameters const & get() const { return m_parameters; } @@ -110,6 +104,16 @@ class LinearSolverParametersInput : public dataRepository::Group /// Adaptive exponent parameter key static constexpr char const * adaptiveExponentString() { return "adaptiveExponent"; } + /// 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 @@ -138,6 +142,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"; } @@ -145,6 +151,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/PhysicsSolverBase.cpp b/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp index 07056cdd12b..834507026a2 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp @@ -42,19 +42,12 @@ PhysicsSolverBase::PhysicsSolverBase( string const & name, m_cflFactor(), m_nextDt( 1e99 ), m_dofManager( name ), - m_usePhysicsScaling(), + m_usePhysicsScaling( 1 ), m_linearSolverParameters( groupKeyStruct::linearSolverParametersString(), this ), m_nonlinearSolverParameters( groupKeyStruct::nonlinearSolverParametersString(), this ), m_solverStatistics( groupKeyStruct::solverStatisticsString(), this ), m_systemSetupTimestamp( 0 ) { - // Physics-scaling is enabled by default only with hypre builds -#ifdef GEOS_USE_HYPRE - integer usePhysicsScaling = 1; -#else - integer usePhysicsScaling = 0; -#endif - setInputFlags( InputFlags::OPTIONAL_NONUNIQUE ); // This sets a flag to indicate that this object is going to select the time step size @@ -101,7 +94,7 @@ PhysicsSolverBase::PhysicsSolverBase( string const & name, setDescription( "Cut time step if linear solution fail without going until max nonlinear iterations." ); registerWrapper( viewKeyStruct::usePhysicsScalingString(), &m_usePhysicsScaling ). - setApplyDefaultValue( usePhysicsScaling ). + setApplyDefaultValue( 1 ). setInputFlag( InputFlags::OPTIONAL ). setDescription( "Enable physics-based scaling of the linear system. Default: true." ); @@ -506,7 +499,7 @@ real64 PhysicsSolverBase::setNextDtBasedOnIterNumber( real64 const & currentDt ) real64 PhysicsSolverBase::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. @@ -552,12 +545,6 @@ real64 PhysicsSolverBase::linearImplicitStep( real64 const & time_n, { Timer timer( m_timers["linear solver total"] ); - // TODO: Trilinos currently requires this, re-evaluate after moving to Tpetra-based solvers - if( m_precond ) - { - m_precond->clear(); - } - { Timer timer_create( m_timers["linear solver create"] ); @@ -566,7 +553,7 @@ real64 PhysicsSolverBase::linearImplicitStep( real64 const & time_n, } // 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 ); @@ -948,7 +935,8 @@ bool PhysicsSolverBase::solveNonlinearSystem( real64 const & time_n, { GEOS_LOG_LEVEL_RANK_0( logInfo::NonlinearSolver, - GEOS_FMT( " Attempt: {:2}, ConfigurationIter: {:2}, NewtonIter: {:2}", dtAttempt, configurationLoopIter, newtonIter )); + GEOS_FMT( " Attempt: {:2}, ConfigurationIter: {:2}, NewtonIter: {:2}", + dtAttempt, configurationLoopIter, newtonIter )); { Timer timer( m_timers["assemble"] ); @@ -1084,12 +1072,6 @@ bool PhysicsSolverBase::solveNonlinearSystem( real64 const & time_n, krylovParams.relTolerance = newtonIter > 0 ? eisenstatWalker( residualNorm, lastResidual, krylovParams ) : krylovParams.weakestTol; } - // TODO: Trilinos currently requires this, re-evaluate after moving to Tpetra-based solvers - if( m_precond ) - { - m_precond->clear(); - } - { Timer timer_setup( m_timers["linear solver create"] ); @@ -1112,7 +1094,9 @@ bool PhysicsSolverBase::solveNonlinearSystem( real64 const & time_n, // Do not allow non converged linear solver - cut time step if( !m_allowNonConvergedLinearSolverSolution && m_linearSolverResult.status == LinearSolverResult::Status::NotConverged ) + { return false; + } } { @@ -1199,6 +1183,16 @@ void PhysicsSolverBase::setupSystem( DomainPartition & domain, solution.create( dofManager.numLocalDofs(), MPI_COMM_GEOS ); } +std::unique_ptr< PreconditionerBase< LAInterface > > +PhysicsSolverBase::createPreconditioner( DomainPartition & GEOS_UNUSED_PARAM( domain ) ) const +{ + // By default, do not create a preconditioner, one will be created internally inside LA backend + return {}; + + // TODO: refactor interfaces to always create preconditioner externally and pass to backends + // return LAInterface::createPreconditioner( m_linearSolverParameters.get() ); +} + void PhysicsSolverBase::assembleSystem( real64 const GEOS_UNUSED_PARAM( time ), real64 const GEOS_UNUSED_PARAM( dt ), DomainPartition & GEOS_UNUSED_PARAM( domain ), @@ -1245,16 +1239,15 @@ 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 ); } if( toFile ) { string const filename = GEOS_FMT( "{}_{:06}_{:02}.mtx", filePrefix.c_str(), cycleNumber, nonlinearIteration ); - obj.write( filename, LAIOutputFormat::NATIVE_ASCII ); - GEOS_LOG_RANK_0( screenName << " written to " << filename ); + obj.write( filename, LAIOutputFormat::MATRIX_MARKET ); + GEOS_LOG_RANK_0( GEOS_FMT( "{} written to {}", screenName, filename ) ); } } @@ -1371,7 +1364,7 @@ void PhysicsSolverBase::solveLinearSystem( DofManager const & dofManager, } GEOS_LOG_LEVEL_RANK_0( logInfo::LinearSolver, - GEOS_FMT( " Last LinSolve(iter,res) = ( {:3}, {:4.2e} )", + GEOS_FMT( " Linear solve: ( iter, res ) = ( {:3}, {:4.2e} )", m_linearSolverResult.numIterations, m_linearSolverResult.residualReduction )); diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp b/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp index 6659092b3d9..7cf0157667b 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp @@ -398,6 +398,14 @@ class PhysicsSolverBase : 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/fluidFlow/SinglePhaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp index 7bc5082406d..a42263183b4 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 "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 "fieldSpecification/LogLevelsInfo.hpp" +#include "mesh/mpiCommunications/CommunicationTools.hpp" #include "physicsSolvers/LogLevelsInfo.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp" @@ -113,6 +113,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 PhysicsSolverBase::createPreconditioner( domain ); + } + } } template< typename BASE > diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp index 2328f0fc403..a85545c0c43 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp @@ -45,6 +45,7 @@ class SinglePhaseFVM : public BASE using BASE::m_discretizationName; 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. @@ -124,6 +125,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/PoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp index ff48cb0b6e5..7db990f9c29 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp @@ -122,7 +122,7 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, // if( !m_precond && m_linearSolverParameters.get().solverType != LinearSolverParameters::SolverType::direct ) // { - // createPreconditioner( domain ); + // m_precond = createPreconditioner( domain ); // } } diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp index cd4b8e796ce..76df6d48faf 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp @@ -59,6 +59,9 @@ SinglePhasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::SinglePhasePoromechan setApplyDefaultValue( 0 ). setInputFlag( InputFlags::OPTIONAL ). setDescription( "The flag to indicate whether a damage solid model is used" ); + + LinearSolverParameters & linearSolverParameters = this->m_linearSolverParameters.get(); + linearSolverParameters.dofsPerNode = 3; } template< typename FLOW_SOLVER, typename MECHANICS_SOLVER > @@ -98,7 +101,7 @@ void SinglePhasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::setupSystem( Dom if( !this->m_precond && this->m_linearSolverParameters.get().solverType != LinearSolverParameters::SolverType::direct ) { - createPreconditioner(); + this->m_precond = createPreconditioner( domain ); } } @@ -119,6 +122,12 @@ void SinglePhasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::initializePostIn getCatalogName(), this->getDataContext(), poromechanicsTargetRegionNames[i], this->flowSolver()->getDataContext() ), InputError ); } + + // Populate sub-block solver parameters for block preconditioner + LinearSolverParameters & linParams = this->m_linearSolverParameters.get(); + linParams.block.resize( 2 ); + linParams.block.subParams[0] = &this->solidMechanicsSolver()->getLinearSolverParameters(); + linParams.block.subParams[1] = &this->flowSolver()->getLinearSolverParameters(); } template<> @@ -379,30 +388,29 @@ void SinglePhasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::assembleElementB } template< typename FLOW_SOLVER, typename MECHANICS_SOLVER > -void SinglePhasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::createPreconditioner() +std::unique_ptr< PreconditionerBase< LAInterface > > +SinglePhasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::createPreconditioner( DomainPartition & domain ) const { - if( this->m_linearSolverParameters.get().preconditionerType == LinearSolverParameters::PreconditionerType::block ) + LinearSolverParameters const & linParams = this->m_linearSolverParameters.get(); + switch( linParams.preconditionerType ) { - auto precond = std::make_unique< BlockPreconditioner< LAInterface > >( BlockShapeOption::UpperTriangular, - SchurComplementOption::RowsumDiagonalProbing, - BlockScalingOption::FrobeniusNorm ); - - auto mechPrecond = LAInterface::createPreconditioner( this->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( this->flowSolver()->getLinearSolverParameters() ); - precond->setupBlock( 1, - { { flow::pressure::key(), { 1, true } } }, - std::move( flowPrecond ) ); + precond->setupBlock( linParams.block.order[toUnderlying( Base::SolverType::SolidMechanics )], + { { solidMechanics::totalDisplacement::key(), { 3, true } } }, + this->solidMechanicsSolver()->createPreconditioner( domain ) ); + precond->setupBlock( linParams.block.order[toUnderlying( Base::SolverType::Flow )], + { { SinglePhaseBase::viewKeyStruct::elemDofFieldString(), { 1, true } } }, + this->flowSolver()->createPreconditioner( domain ) ); - this->m_precond = std::move( precond ); - } - else - { - //TODO: Revisit this part such that is coherent across physics solver - //m_precond = LAInterface::createPreconditioner( m_linearSolverParameters.get() ); + return precond; + } + default: + { + return PhysicsSolverBase::createPreconditioner( domain ); + } } } diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp index 66d12a7da9c..c0232c8ef25 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp @@ -51,9 +51,6 @@ class SinglePhasePoromechanics : public PoromechanicsSolver< FLOW_SOLVER, MECHAN SinglePhasePoromechanics( const string & name, dataRepository::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. @@ -94,6 +91,9 @@ class SinglePhasePoromechanics : public PoromechanicsSolver< FLOW_SOLVER, MECHAN 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, @@ -133,9 +133,6 @@ class SinglePhasePoromechanics : public PoromechanicsSolver< FLOW_SOLVER, MECHAN * @param[in] subRegion the element subRegion */ virtual void updateBulkDensity( ElementSubRegionBase & subRegion ) override; - - void createPreconditioner(); - }; } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index 0f594f2703b..c8654d14c52 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -26,6 +26,7 @@ #include "kernels/ExplicitFiniteStrain.hpp" #include "kernels/FixedStressThermoPoromechanics.hpp" +#include "common/GEOS_RAJA_Interface.hpp" #include "constitutive/ConstitutiveManager.hpp" #include "common/GEOS_RAJA_Interface.hpp" #include "discretizationMethods/NumericalMethodsManager.hpp" @@ -33,7 +34,9 @@ #include "fieldSpecification/TractionBoundaryCondition.hpp" #include "finiteElement/FiniteElementDiscretizationManager.hpp" #include "finiteElement/FiniteElementDiscretization.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" @@ -129,27 +132,21 @@ SolidMechanicsLagrangianFEM::SolidMechanicsLagrangianFEM( const string & name, setInputFlag( InputFlags::FALSE ). setDescription( "The maximum force contribution in the problem domain." ); -} - -void SolidMechanicsLagrangianFEM::postInputInitialization() -{ - PhysicsSolverBase::postInputInitialization(); - - // configure AMG + // Set physics-dependent parameters for linear solver LinearSolverParameters & linParams = m_linearSolverParameters.get(); - linParams.isSymmetric = true; linParams.dofsPerNode = 3; linParams.amg.numFunctions = linParams.dofsPerNode; - linParams.amg.separateComponents = true; - m_surfaceGenerator = this->getParent().getGroupPointer< PhysicsSolverBase >( m_surfaceGeneratorName ); + // Must change default value to avoid being overwritten if attribute not specified in XML + m_linearSolverParameters.getWrapper< integer >( LinearSolverParametersInput::viewKeyStruct::amgSeparateComponentsString() ).setApplyDefaultValue( true ); } -SolidMechanicsLagrangianFEM::~SolidMechanicsLagrangianFEM() +void SolidMechanicsLagrangianFEM::postInputInitialization() { - // TODO Auto-generated destructor stub -} + PhysicsSolverBase::postInputInitialization(); + m_surfaceGenerator = this->getParent().getGroupPointer< PhysicsSolverBase >( m_surfaceGeneratorName ); +} void SolidMechanicsLagrangianFEM::registerDataOnMesh( Group & meshBodies ) { @@ -1010,6 +1007,47 @@ 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 ); + } +} + +void SolidMechanicsLagrangianFEM::computeRigidBodyModes( DomainPartition & domain ) const +{ + LinearSolverParameters const & linParams = m_linearSolverParameters.get(); + if( linParams.amg.nullSpaceType == LinearSolverParameters::AMG::NullSpaceType::rigidBodyModes && m_rigidBodyModes.empty() ) + { + 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() ) ); + } + } ); + } +} + +std::unique_ptr< PreconditionerBase< LAInterface > > +SolidMechanicsLagrangianFEM::createPreconditioner( DomainPartition & domain ) const +{ + LinearSolverParameters const & linParams = m_linearSolverParameters.get(); + + switch( linParams.preconditionerType ) + { + default: + { + return PhysicsSolverBase::createPreconditioner( domain ); + } + } } void SolidMechanicsLagrangianFEM::assembleSystem( real64 const GEOS_UNUSED_PARAM( time_n ), @@ -1249,7 +1287,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 68c52e3a82f..89579b79d15 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp @@ -65,18 +65,6 @@ class SolidMechanicsLagrangianFEM : public PhysicsSolverBase 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 PhysicsSolverBase::CatalogInterface::CatalogType */ @@ -125,6 +113,9 @@ class SolidMechanicsLagrangianFEM : public PhysicsSolverBase 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, @@ -265,13 +256,11 @@ class SolidMechanicsLagrangianFEM : public PhysicsSolverBase real64 & getMaxForce() { return m_maxForce; } real64 const & getMaxForce() const { return m_maxForce; } - arrayView1d< ParallelVector > const & getRigidBodyModes() const - { - return m_rigidBodyModes; - } + void computeRigidBodyModes( DomainPartition & domain ) const; - array1d< ParallelVector > & getRigidBodyModes() + arrayView1d< ParallelVector > const & getRigidBodyModes( DomainPartition & domain ) const { + computeRigidBodyModes( domain ); return m_rigidBodyModes; } @@ -307,8 +296,8 @@ class SolidMechanicsLagrangianFEM : public PhysicsSolverBase /// Flag to indicate that the solver is going to perform stress initialization bool m_performStressInitialization; - /// Rigid body modes - array1d< ParallelVector > m_rigidBodyModes; + /// Rigid body modes; TODO remove mutable hack + mutable array1d< ParallelVector > m_rigidBodyModes; real64 m_contactPenaltyStiffness; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp index 8e1c3a159eb..2911d338f9a 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp @@ -872,25 +872,28 @@ real64 SolidMechanicsLagrangeContact::calculateContactResidualNorm( DomainPartit return sqrt( stickResidual * stickResidual + slipResidual * slipResidual + openResidual * openResidual ); } -void SolidMechanicsLagrangeContact::createPreconditioner( DomainPartition const & domain ) +std::unique_ptr< PreconditionerBase< LAInterface > > +SolidMechanicsLagrangeContact::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 = 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 > >(); @@ -901,9 +904,8 @@ void SolidMechanicsLagrangeContact::createPreconditioner( DomainPartition const } 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 @@ -916,31 +918,20 @@ void SolidMechanicsLagrangeContact::createPreconditioner( DomainPartition const { { contact::traction::key(), { 3, true } } }, std::move( tracPrecond ) ); - if( mechParams.amg.nullSpaceType == LinearSolverParameters::AMG::NullSpaceType::rigidBodyModes ) - { - if( getRigidBodyModes().empty() ) - { - MeshLevel const & mesh = domain.getMeshBody( 0 ).getBaseDiscretization(); - LAIHelperFunctions::computeRigidBodyModes( mesh, - m_dofManager, - { solidMechanics::totalDisplacement::key() }, - getRigidBodyModes() ); - } - } - // Preconditioner for the Schur complement: mechPrecond - std::unique_ptr< PreconditionerBase< LAInterface > > mechPrecond = LAInterface::createPreconditioner( mechParams, getRigidBodyModes() ); + std::unique_ptr< PreconditionerBase< LAInterface > > mechPrecond = LAInterface::createPreconditioner( mechParams, getRigidBodyModes( domain ) ); precond->setupBlock( 1, - { { solidMechanics::totalDisplacement::key(), { 3, true } } }, + { { fields::solidMechanics::totalDisplacement::key(), { 3, true } } }, std::move( mechPrecond ) ); - 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 SolidMechanicsLagrangeContact::computeRotationMatrices( DomainPartition & domain ) const diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp index 2f7104dbc52..06b6bf97d36 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp @@ -65,6 +65,9 @@ class SolidMechanicsLagrangeContact : public ContactSolverBase ParallelVector & solution, bool const setSparsity = true ) override final; + virtual std::unique_ptr< PreconditionerBase< LAInterface > > + createPreconditioner( DomainPartition & domain ) const override; + virtual void implicitStepSetup( real64 const & time_n, real64 const & dt, @@ -187,8 +190,6 @@ class SolidMechanicsLagrangeContact : public ContactSolverBase static const localIndex m_maxFaceNodes; // Maximum number of nodes on a contact face - void createPreconditioner( DomainPartition const & domain ); - void computeFaceDisplacementJump( DomainPartition & domain ); struct viewKeyStruct : ContactSolverBase::viewKeyStruct diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index 81dbf03b5fc..875611d0423 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -1973,6 +1973,9 @@ Information output from lower logLevels is added with the desired log level + + + @@ -1983,7 +1986,7 @@ Information output from lower logLevels is added with the desired log level - + @@ -1991,8 +1994,12 @@ Information output from lower logLevels is added with the desired log level - + + + + + @@ -2001,13 +2008,17 @@ Information output from lower logLevels is added with the desired log level - + + + + + - + @@ -2017,7 +2028,7 @@ Information output from lower logLevels is added with the desired log level - + @@ -2044,13 +2055,38 @@ Information output from lower logLevels is added with the desired log level - Linear solver information - Print linear solver configuration--> - + - + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -2078,7 +2114,7 @@ Information output from lower logLevels is added with the desired log level - + @@ -2093,12 +2129,12 @@ Information output from lower logLevels is added with the desired log level - + - + diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other index b8d0088321c..69761fa702a 100644 --- a/src/coreComponents/schema/schema.xsd.other +++ b/src/coreComponents/schema/schema.xsd.other @@ -451,7 +451,12 @@ - + + + + + + diff --git a/src/coreComponents/unitTests/linearAlgebraTests/testDofManager.cpp b/src/coreComponents/unitTests/linearAlgebraTests/testDofManager.cpp index c7a0dc9212d..dd6de25cdf3 100644 --- a/src/coreComponents/unitTests/linearAlgebraTests/testDofManager.cpp +++ b/src/coreComponents/unitTests/linearAlgebraTests/testDofManager.cpp @@ -589,7 +589,6 @@ void DofManagerSparsityTest< LAI >::test( stdVector< FieldDesc > fields, CRSMatrix< real64, globalIndex > localMatrix; localMatrix.assimilate< parallelHostPolicy >( std::move( localPattern ) ); pattern.create( localMatrix.toViewConst(), dofManager.numLocalDofs(), MPI_COMM_GEOS ); - pattern.set( 1.0 ); } CRSMatrix< real64, globalIndex > localPatternExpected( numLocalDof, @@ -599,7 +598,6 @@ void DofManagerSparsityTest< LAI >::test( stdVector< FieldDesc > fields, for( FieldDesc const & f : fields ) { - GEOS_LOG_RANK( "rankOffset = "<