Skip to content

GitLab

Commit 1b84cfb5 authored by Nils Kohl's avatar Nils Kohl 🌝
Browse files

Merge branch 'kohl/petsc-dynamic-sparse-matrix-allocation' into 'master' 

PETSc dynamic sparse matrix allocation

See merge request !476
parents cad4b9db 83b8b855
Pipeline #37059 passed with stages
in 147 minutes and 42 seconds
Show whitespace changes
Inline Side-by-side
apps/annulus.cpp
View file @ 1b84cfb5
......@@ -91,7 +91,7 @@ void solveProblem( std::shared_ptr< hyteg::PrimitiveStorage >& storage, uint_t l
uint_t localDoFs = numberOfLocalDoFs< enumTag >( *storage, level );
// assemble matrices
PETScSparseMatrix< opType > lapMat( localDoFs, globalDoFs );
PETScSparseMatrix< opType > lapMat( "Mat" );
switch ( verbosity )
{
......
apps/benchmarks/PetscCompare-2D-P2-Apply/PetscCompare-2D-P2-Apply.cpp
View file @ 1b84cfb5
......@@ -126,18 +126,17 @@ int main( int argc, char* argv[] )
wcTimingTreeApp.stop( "Interpolation" );
wcTimingTreeApp.start( "Enumeration" );
const uint_t globalDoFs = hyteg::numberOfGlobalDoFs< P2FunctionTag >( *storage, level );
const uint_t localDoFs = hyteg::numberOfLocalDoFs< P2FunctionTag >( *storage, level );
numerator.enumerate( level );
wcTimingTreeApp.stop( "Enumeration" );
LIKWID_MARKER_START( "PETSc-setup" );
wcTimingTreeApp.start( "Petsc setup" );
hyteg::PETScSparseMatrix< hyteg::P2ConstantLaplaceOperator > matPetsc( localDoFs, globalDoFs );
hyteg::PETScSparseMatrix< hyteg::P2ConstantLaplaceOperator > matPetsc;
matPetsc.createMatrixFromOperator( mass, level, numerator, hyteg::Inner );
hyteg::PETScVector< real_t, hyteg::P2Function > vecPetsc( localDoFs );
hyteg::PETScVector< real_t, hyteg::P2Function > vecPetsc;
vecPetsc.createVectorFromFunction( x, numerator, level, hyteg::Inner );
hyteg::PETScVector< real_t, hyteg::P2Function > dstvecPetsc( localDoFs );
hyteg::PETScVector< real_t, hyteg::P2Function > dstvecPetsc;
dstvecPetsc.createVectorFromFunction( x, numerator, level, hyteg::Inner );
wcTimingTreeApp.stop( "Petsc setup" );
LIKWID_MARKER_STOP( "PETSc-setup" );
......
apps/benchmarks/PetscCompare-3D-P1-Apply/PetscCompare-3D-P1-Apply.cpp
View file @ 1b84cfb5
......@@ -125,18 +125,17 @@ int main( int argc, char* argv[] )
wcTimingTreeApp.stop( "Interpolation" );
wcTimingTreeApp.start( "Enumeration" );
const uint_t globalDoFs = hyteg::numberOfGlobalDoFs< P1FunctionTag >( *storage, level );
const uint_t localDoFs = hyteg::numberOfLocalDoFs< P1FunctionTag >( *storage, level );
numerator.enumerate( level );
wcTimingTreeApp.stop( "Enumeration" );
LIKWID_MARKER_START( "PETSc-setup" );
wcTimingTreeApp.start( "Petsc setup" );
hyteg::PETScSparseMatrix< hyteg::P1ConstantLaplaceOperator > matPetsc( localDoFs, globalDoFs );
hyteg::PETScSparseMatrix< hyteg::P1ConstantLaplaceOperator > matPetsc;
matPetsc.createMatrixFromOperator( mass, level, numerator, hyteg::Inner );
hyteg::PETScVector< real_t, hyteg::P1Function > vecPetsc( localDoFs );
hyteg::PETScVector< real_t, hyteg::P1Function > vecPetsc;
vecPetsc.createVectorFromFunction( x, numerator, level, hyteg::Inner );
hyteg::PETScVector< real_t, hyteg::P1Function > dstvecPetsc( localDoFs );
hyteg::PETScVector< real_t, hyteg::P1Function > dstvecPetsc;
dstvecPetsc.createVectorFromFunction( x, numerator, level, hyteg::Inner );
wcTimingTreeApp.stop( "Petsc setup" );
LIKWID_MARKER_STOP( "PETSc-setup" );
......
apps/benchmarks/PetscCompare/PetscCompare.cpp
View file @ 1b84cfb5
......@@ -56,12 +56,13 @@ int main( int argc, char* argv[] )
PETScManager petscManager( &argc, &argv );
auto cfg = std::make_shared< walberla::config::Config >();
if( env.config() == nullptr )
if ( env.config() == nullptr )
{
auto defaultFile = "./PetscCompare.prm";
WALBERLA_LOG_PROGRESS_ON_ROOT( "No Parameter file given loading default parameter file: " << defaultFile );
cfg->readParameterFile( defaultFile );
} else
}
else
{
cfg = env.config();
}
......@@ -69,14 +70,15 @@ int main( int argc, char* argv[] )
const uint_t level = mainConf.getParameter< uint_t >( "level" );
std::shared_ptr< hyteg::MeshInfo > meshInfo;
if( mainConf.getParameter< bool >( "useMeshFile" ) )
if ( mainConf.getParameter< bool >( "useMeshFile" ) )
{
std::string meshFileName = mainConf.getParameter< std::string >( "mesh" );
meshInfo = std::make_shared< hyteg::MeshInfo >( hyteg::MeshInfo::fromGmshFile( meshFileName ) );
} else
}
else
{
uint_t numberOfFaces = mainConf.getParameter< uint_t >( "numberOfFaces" );
if( mainConf.getParameter< bool >( "facesTimesProcs" ) )
if ( mainConf.getParameter< bool >( "facesTimesProcs" ) )
{
meshInfo = std::make_shared< hyteg::MeshInfo >(
hyteg::MeshInfo::meshFaceChain( numberOfFaces * uint_c( walberla::MPIManager::instance()->numProcesses() ) ) );
......@@ -101,7 +103,7 @@ int main( int argc, char* argv[] )
std::shared_ptr< hyteg::PrimitiveStorage > storage = std::make_shared< hyteg::PrimitiveStorage >( setupStorage, timingTree );
if( mainConf.getParameter< bool >( "writeDomain" ) )
if ( mainConf.getParameter< bool >( "writeDomain" ) )
{
hyteg::writeDomainPartitioningVTK( storage, "./output", "domain" );
}
......@@ -162,11 +164,12 @@ int main( int argc, char* argv[] )
numerator.enumerate( level );
LIKWID_MARKER_START( "PETSc-setup" );
hyteg::PETScSparseMatrix< hyteg::P1ConstantLaplaceOperator > matPetsc( localDoFs, totalDoFs );
hyteg::PETScSparseMatrix< hyteg::P1ConstantLaplaceOperator > matPetsc;
matPetsc.createMatrixFromOperator( mass, level, numerator, hyteg::Inner );
hyteg::PETScVector< real_t, hyteg::P1Function > vecPetsc( localDoFs );
hyteg::PETScVector< real_t, hyteg::P1Function > vecPetsc;
vecPetsc.createVectorFromFunction( x, numerator, level, hyteg::Inner );
hyteg::PETScVector< real_t, hyteg::P1Function > dstvecPetsc( localDoFs );
hyteg::PETScVector< real_t, hyteg::P1Function > dstvecPetsc;
dstvecPetsc.createVectorFromFunction( x, numerator, level, hyteg::Inner );
LIKWID_MARKER_STOP( "PETSc-setup" );
walberla::WcTimer timer;
......@@ -196,12 +199,12 @@ int main( int argc, char* argv[] )
//dstvecPetsc.print("../output/vector1.vec");
diff.assign( {1.0, -1.0}, {z, y}, level, hyteg::All );
diff.assign( { 1.0, -1.0 }, { z, y }, level, hyteg::All );
if( mainConf.getParameter< bool >( "VTKOutput" ) )
if ( mainConf.getParameter< bool >( "VTKOutput" ) )
{
WALBERLA_LOG_INFO_ON_ROOT( "writing VTK output" );
hyteg::VTKOutput vtkOutput("./output", "petscCompare", storage);
hyteg::VTKOutput vtkOutput( "./output", "petscCompare", storage );
vtkOutput.add( x );
vtkOutput.add( z );
vtkOutput.add( y );
......@@ -212,12 +215,12 @@ int main( int argc, char* argv[] )
walberla::WcTimingTree tt = timingTree->getReduced();
auto tt2 = tt.getCopyWithRemainder();
if( mainConf.getParameter< bool >( "printTiming" ) )
if ( mainConf.getParameter< bool >( "printTiming" ) )
{
WALBERLA_LOG_INFO_ON_ROOT( tt2 );
}
if( mainConf.getParameter< bool >( "writeJSON" ) )
if ( mainConf.getParameter< bool >( "writeJSON" ) )
{
nlohmann::json ttjson = nlohmann::json( tt2 );
std::ofstream o( "PetscCompareOutput.json" );
......
src/hyteg/petsc/PETScBlockPreconditionedStokesSolver.hpp
View file @ 1b84cfb5
/*
* Copyright (c) 2017-2022 Nils Kohl.
*
* This file is part of HyTeG
* (see https://i10git.cs.fau.de/hyteg/hyteg).
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#pragma once
#include <memory>
......@@ -39,21 +58,12 @@ class PETScBlockPreconditionedStokesSolver : public Solver< OperatorType >
: allocatedLevel_( level )
, petscCommunicator_( storage->getSplitCommunicatorByPrimitiveDistribution() )
, num( "numerator", storage, level, level )
, Amat( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ),
numberOfGlobalDoFs< typename FunctionType::Tag >( *storage, level, petscCommunicator_ ),
"Amat",
petscCommunicator_ )
, AmatNonEliminatedBC( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ),
numberOfGlobalDoFs< typename FunctionType::Tag >( *storage, level, petscCommunicator_ ),
"AmatNonEliminatedBC",
petscCommunicator_ )
, Pmat( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ),
numberOfGlobalDoFs< typename FunctionType::Tag >( *storage, level, petscCommunicator_ ),
"Pmat",
petscCommunicator_ )
, xVec( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ), "xVec", petscCommunicator_ )
, bVec( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ), "bVec", petscCommunicator_ )
, nullspaceVec_( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ), "nullspaceVec", petscCommunicator_ )
, Amat( "Amat", petscCommunicator_ )
, AmatNonEliminatedBC( "AmatNonEliminatedBC", petscCommunicator_ )
, Pmat( "Pmat", petscCommunicator_ )
, xVec( "xVec", petscCommunicator_ )
, bVec( "bVec", petscCommunicator_ )
, nullspaceVec_( "nullspaceVec", petscCommunicator_ )
, storage_( storage )
, tolerance_( tolerance )
, maxIterations_( maxIterations )
......@@ -208,8 +218,7 @@ class PETScBlockPreconditionedStokesSolver : public Solver< OperatorType >
case 1:
PCSetType( pc_u, PCJACOBI );
break;
case 3:
{
case 3: {
PCSetType( pc_u, PCHYPRE );
break;
}
......
src/hyteg/petsc/PETScCGSolver.hpp
View file @ 1b84cfb5
/*
* Copyright (c) 2017-2019 Dominik Thoennes, Nils Kohl.
* Copyright (c) 2017-2022 Dominik Thoennes, Nils Kohl.
*
* This file is part of HyTeG
* (see https://i10git.cs.fau.de/hyteg/hyteg).
......@@ -44,17 +44,11 @@ class PETScCGSolver : public Solver< OperatorType >
: allocatedLevel_( level )
, petscCommunicator_( storage->getSplitCommunicatorByPrimitiveDistribution() )
, num( "numerator", storage, level, level )
, Amat( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ),
numberOfGlobalDoFs< typename FunctionType::Tag >( *storage, level, petscCommunicator_ ),
"Amat",
petscCommunicator_ )
, AmatNonEliminatedBC( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ),
numberOfGlobalDoFs< typename FunctionType::Tag >( *storage, level, petscCommunicator_ ),
"AmatNonEliminatedBC",
petscCommunicator_ )
, xVec( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ), "xVec", petscCommunicator_ )
, bVec( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ), "bVec", petscCommunicator_ )
, nullspaceVec_( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ), "nullspaceVec", petscCommunicator_ )
, Amat( "Amat", petscCommunicator_ )
, AmatNonEliminatedBC( "AmatNonEliminatedBC", petscCommunicator_ )
, xVec( "xVec", petscCommunicator_ )
, bVec( "bVec", petscCommunicator_ )
, nullspaceVec_( "nullspaceVec", petscCommunicator_ )
, flag_( hyteg::All )
, nullSpaceSet_( false )
, reassembleMatrix_( false )
......
src/hyteg/petsc/PETScExportLinearSystem.hpp
View file @ 1b84cfb5
/*
* Copyright (c) 2017-2019 Nils Kohl, Dominik Thoennes, Marcus Mohr.
* Copyright (c) 2017-2022 Nils Kohl, Dominik Thoennes, Marcus Mohr.
*
* This file is part of HyTeG
* (see https://i10git.cs.fau.de/hyteg/hyteg).
......@@ -92,7 +92,7 @@ void exportLinearSystem( OperatorType op,
{
WALBERLA_LOG_INFO_ON_ROOT( " * Converting Operator to PETSc matrix" );
}
PETScSparseMatrix< OperatorType > petscMatrix( localDoFs, globalDoFs, nameMatrix.c_str() );
PETScSparseMatrix< OperatorType > petscMatrix( nameMatrix.c_str() );
FunctionType< idx_t > numerator( "numerator", storage, level, level );
numerator.enumerate( level );
petscMatrix.createMatrixFromOperator( op, level, numerator );
......
src/hyteg/petsc/PETScExportOperatorMatrix.hpp
View file @ 1b84cfb5
......@@ -115,7 +115,7 @@ void exportOperator( OperatorType& op,
{
WALBERLA_LOG_INFO_ON_ROOT( " * Converting Operator to PETSc matrix" )
}
PETScSparseMatrix< OperatorType > petscMatrix( localDoFs, globalDoFs, matrixName.c_str() );
PETScSparseMatrix< OperatorType > petscMatrix( matrixName.c_str() );
typename OperatorType::srcType::template FunctionType< idx_t > numerator( "numerator", storage, level, level );
numerator.enumerate( level );
petscMatrix.createMatrixFromOperator( op, level, numerator );
......
src/hyteg/petsc/PETScLUSolver.hpp
View file @ 1b84cfb5
/*
* Copyright (c) 2017-2019 Boerge Struempfel, Daniel Drzisga, Dominik Thoennes, Nils Kohl.
* Copyright (c) 2017-2022 Boerge Struempfel, Daniel Drzisga, Dominik Thoennes, Nils Kohl.
*
* This file is part of HyTeG
* (see https://i10git.cs.fau.de/hyteg/hyteg).
......@@ -53,20 +53,11 @@ class PETScLUSolver : public Solver< OperatorType >
, allocatedLevel_( level )
, petscCommunicator_( storage->getSplitCommunicatorByPrimitiveDistribution() )
, num( "numerator", storage, level, level )
, Amat( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ),
numberOfGlobalDoFs< typename FunctionType::Tag >( *storage, level, petscCommunicator_ ),
"Amat",
petscCommunicator_ )
, AmatUnsymmetric( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ),
numberOfGlobalDoFs< typename FunctionType::Tag >( *storage, level, petscCommunicator_ ),
"AmatUnsymmetric",
petscCommunicator_ )
, AmatTmp( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ),
numberOfGlobalDoFs< typename FunctionType::Tag >( *storage, level, petscCommunicator_ ),
"AmatTmp",
petscCommunicator_ )
, xVec( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ), "xVec", petscCommunicator_ )
, bVec( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ), "bVec", petscCommunicator_ )
, Amat( "Amat", petscCommunicator_ )
, AmatUnsymmetric( "AmatUnsymmetric", petscCommunicator_ )
, AmatTmp( "AmatTmp", petscCommunicator_ )
, xVec( "xVec", petscCommunicator_ )
, bVec( "bVec", petscCommunicator_ )
#if 0
, inKernel( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ) )
#endif
......@@ -231,8 +222,9 @@ class PETScLUSolver : public Solver< OperatorType >
storage_->getTimingTree()->start( "RHS vector setup" );
b.assign( {1.0}, {x}, level, DirichletBoundary );
b.assign( { 1.0 }, { x }, level, DirichletBoundary );
bVec.createVectorFromFunction( b, num, level, All );
xVec.createVectorFromFunction( x, num, level, All );
if ( assumeSymmetry_ )
{
......
src/hyteg/petsc/PETScMinResSolver.hpp
View file @ 1b84cfb5
/*
* Copyright (c) 2017-2019 Dominik Thoennes, Nils Kohl.
* Copyright (c) 2017-2022 Dominik Thoennes, Nils Kohl.
*
* This file is part of HyTeG
* (see https://i10git.cs.fau.de/hyteg/hyteg).
......@@ -44,17 +44,11 @@ class PETScMinResSolver : public Solver< OperatorType >
: allocatedLevel_( level )
, petscCommunicator_( storage->getSplitCommunicatorByPrimitiveDistribution() )
, num( "numerator", storage, level, level )
, Amat( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ),
numberOfGlobalDoFs< typename FunctionType::Tag >( *storage, level, petscCommunicator_ ),
"Amat",
petscCommunicator_ )
, AmatNonEliminatedBC( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ),
numberOfGlobalDoFs< typename FunctionType::Tag >( *storage, level, petscCommunicator_ ),
"AmatNonEliminatedBC",
petscCommunicator_ )
, xVec( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ), "xVec", petscCommunicator_ )
, bVec( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ), "bVec", petscCommunicator_ )
, nullspaceVec_( numberOfLocalDoFs< typename FunctionType::Tag >( *storage, level ), "nullspaceVec", petscCommunicator_ )
, Amat( "Amat", petscCommunicator_ )
, AmatNonEliminatedBC( "AmatNonEliminatedBC", petscCommunicator_ )
, xVec( "xVec", petscCommunicator_ )
, bVec( "bVec", petscCommunicator_ )
, nullspaceVec_( "nullspaceVec", petscCommunicator_ )
, flag_( hyteg::All )
, nullSpaceSet_( false )
, reassembleMatrix_( false )
......
src/hyteg/petsc/PETScSparseMatrix.hpp
View file @ 1b84cfb5
/*
* Copyright (c) 2017-2019 Boerge Struempfel, Daniel Drzisga, Dominik Thoennes, Marcus Mohr, Nils Kohl.
* Copyright (c) 2017-2022 Boerge Struempfel, Daniel Drzisga, Dominik Thoennes, Marcus Mohr, Nils Kohl.
*
* This file is part of HyTeG
* (see https://i10git.cs.fau.de/hyteg/hyteg).
......@@ -51,78 +51,34 @@ class PETScSparseMatrix
template < typename ValueType >
using FunctionTypeDst = typename OperatorType::dstType::template FunctionType< ValueType >;
PETScSparseMatrix() = delete;
PETScSparseMatrix( uint_t localRows,
uint_t localCols,
uint_t globalRows,
uint_t globalCols,
const char name[] = "Mat",
PETScSparseMatrix( const std::string name = "Mat",
const MPI_Comm& petscCommunicator = walberla::mpi::MPIManager::instance()->comm() )
: petscCommunicator_( petscCommunicator )
: name_( name )
, petscCommunicator_( petscCommunicator )
, allocated_( false )
, assembled_( false )
{
MatCreate( petscCommunicator, &mat );
MatSetType( mat, MATMPIAIJ );
MatSetSizes( mat,
static_cast< PetscInt >( localRows ),
static_cast< PetscInt >( localCols ),
static_cast< PetscInt >( globalRows ),
static_cast< PetscInt >( globalCols ) );
// Roughly overestimate number of non-zero entries for faster assembly of matrix
MatMPIAIJSetPreallocation( mat, 500, NULL, 500, NULL );
setName( name );
reset();
}
PETScSparseMatrix( uint_t localSize,
uint_t globalSize,
const char name[] = "Mat",
const MPI_Comm& petscCommunicator = walberla::mpi::MPIManager::instance()->comm() )
: PETScSparseMatrix( localSize, localSize, globalSize, globalSize, name, petscCommunicator )
{}
PETScSparseMatrix( const std::shared_ptr< PrimitiveStorage >& storage,
const uint_t& level,
const char name[] = "Mat",
const MPI_Comm& petscCommunicator = walberla::mpi::MPIManager::instance()->comm() )
: PETScSparseMatrix( numberOfLocalDoFs< typename OperatorType::dstType::Tag >( *storage, level ),
numberOfGlobalDoFs< typename OperatorType::dstType::Tag >( *storage, level, petscCommunicator ),
name,
petscCommunicator )
{}
PETScSparseMatrix( const FunctionTypeSrc< idx_t >& enumerator,
const uint_t& level,
const char name[] = "Mat",
const MPI_Comm& petscCommunicator = walberla::mpi::MPIManager::instance()->comm() )
: PETScSparseMatrix( numberOfLocalDoFs( enumerator, level ),
numberOfGlobalDoFs( enumerator, level, petscCommunicator ),
name,
petscCommunicator )
{}
PETScSparseMatrix( const FunctionTypeSrc< idx_t >& enumeratorSrc,
const FunctionTypeDst< idx_t >& enumeratorDst,
const uint_t& level,
const char name[] = "Mat",
const MPI_Comm& petscCommunicator = walberla::mpi::MPIManager::instance()->comm() )
: PETScSparseMatrix( numberOfLocalDoFs( enumeratorDst, level ),
numberOfLocalDoFs( enumeratorSrc, level ),
numberOfGlobalDoFs( enumeratorDst, level, petscCommunicator ),
numberOfGlobalDoFs( enumeratorSrc, level, petscCommunicator ),
name,
petscCommunicator )
{}
virtual ~PETScSparseMatrix() { MatDestroy( &mat ); }
virtual ~PETScSparseMatrix()
{
if ( allocated_ )
{
MatDestroy( &mat );
}
}
inline void createMatrixFromOperator( const OperatorType& op,
uint_t level,
const FunctionTypeSrc< idx_t >& numerator,
DoFType flag = All )
{
const uint_t localRows = numberOfLocalDoFs( numerator, level );
const uint_t localCols = numberOfLocalDoFs( numerator, level );
const uint_t globalRows = numberOfGlobalDoFs( numerator, level, petscCommunicator_ );
const uint_t globalCols = numberOfGlobalDoFs( numerator, level, petscCommunicator_ );
allocateSparseMatrix( localRows, localCols, globalRows, globalCols );
auto proxy = std::make_shared< PETScSparseMatrixProxy >( mat );
op.toMatrix( proxy, numerator, numerator, level, flag );
......@@ -137,6 +93,13 @@ class PETScSparseMatrix
const FunctionTypeDst< idx_t >& numeratorDst,
DoFType flag = All )
{
const uint_t localRows = numberOfLocalDoFs( numeratorDst, level );
const uint_t localCols = numberOfLocalDoFs( numeratorSrc, level );
const uint_t globalRows = numberOfGlobalDoFs( numeratorDst, level, petscCommunicator_ );
const uint_t globalCols = numberOfGlobalDoFs( numeratorSrc, level, petscCommunicator_ );
allocateSparseMatrix( localRows, localCols, globalRows, globalCols );
auto proxy = std::make_shared< PETScSparseMatrixProxy >( mat );
op.toMatrix( proxy, numeratorSrc, numeratorDst, level, flag );
......@@ -263,7 +226,13 @@ class PETScSparseMatrix
inline void reset() { assembled_ = false; }
/// \brief Sets all entries of the matrix to zero.
inline void zeroEntries() { MatZeroEntries( mat ); }
inline void zeroEntries()
{
if ( allocated_ )
{
MatZeroEntries( mat );
}
}
inline void setName( const char name[] ) { PetscObjectSetName( (PetscObject) mat, name ); }
......@@ -322,9 +291,32 @@ class PETScSparseMatrix
};
protected:
std::string name_;
MPI_Comm petscCommunicator_;
Mat mat;
bool allocated_;
bool assembled_;
private:
inline void allocateSparseMatrix( uint_t localRows, uint_t localCols, uint_t globalRows, uint_t globalCols )
{
if ( !allocated_ )
{
MatCreate( petscCommunicator_, &mat );
MatSetType( mat, MATMPIAIJ );
MatSetSizes( mat,
static_cast< PetscInt >( localRows ),
static_cast< PetscInt >( localCols ),
static_cast< PetscInt >( globalRows ),
static_cast< PetscInt >( globalCols ) );
// Roughly overestimate number of non-zero entries for faster assembly of matrix
MatMPIAIJSetPreallocation( mat, 500, NULL, 500, NULL );
setName( name_.c_str() );
reset();
allocated_ = true;
}
}
};
} // namespace hyteg
......
src/hyteg/petsc/PETScVector.hpp
View file @ 1b84cfb5
/*
* Copyright (c) 2017-2019 Boerge Struempfel, Daniel Drzisga, Dominik Thoennes, Nils Kohl.
* Copyright (c) 2017-2022 Boerge Struempfel, Daniel Drzisga, Dominik Thoennes, Nils Kohl.
*
* This file is part of HyTeG
* (see https://i10git.cs.fau.de/hyteg/hyteg).
......@@ -41,7 +41,12 @@ template < typename ValueType, template < class > class FunctionType >
class PETScVector
{
public:
PETScVector() = delete;
PETScVector( const std::string& name = "Vec",
const MPI_Comm& petscCommunicator = walberla::mpi::MPIManager::instance()->comm() )
: name_( name )
, petscCommunicator_( petscCommunicator )
, allocated_( false )
{}
PETScVector( const FunctionType< ValueType >& function,
const FunctionType< idx_t >& numerator,
......@@ -49,24 +54,18 @@ class PETScVector
const DoFType& flag = All,
const std::string& name = "Vec",
const MPI_Comm& petscCommunicator = walberla::mpi::MPIManager::instance()->comm() )
: PETScVector( numberOfLocalDoFs( function, level ), name, petscCommunicator )
: PETScVector( name, petscCommunicator )
{
createVectorFromFunction( function, numerator, level, flag );
}
PETScVector( uint_t localSize,
const std::string& name = "Vec",
const MPI_Comm& petscCommunicator = walberla::mpi::MPIManager::instance()->comm() )
: petscCommunicator_( petscCommunicator )
~PETScVector()
{
VecCreate( petscCommunicator, &vec );
VecSetType( vec, VECSTANDARD );
VecSetSizes( vec, (idx_t) localSize, PETSC_DECIDE );
VecSetUp( vec );
setName( name.c_str() );
if ( allocated_ )
{
VecDestroy( &vec );
}
}
~PETScVector() { VecDestroy( &vec ); }
MPI_Comm getCommunicator() const { return petscCommunicator_; }
......@@ -75,6 +74,9 @@ class PETScVector
uint_t level,
DoFType flag = All )
{
const auto localSize = numberOfLocalDoFs( src, level );
allocateVector( localSize );
auto proxy = std::make_shared< PETScVectorProxy >( vec );
hyteg::createVectorFromFunction( src, numerator, proxy, level, flag );
......@@ -87,6 +89,8 @@ class PETScVector
uint_t level,
DoFType flag = All )
{
WALBERLA_CHECK( allocated_, "Cannot create function from PETSc vector - the vector wasn't even allocated." );
auto proxy = std::make_shared< PETScVectorProxy >( vec );
hyteg::createFunctionFromVector( dst, numerator, proxy, level, flag );
}
......@@ -112,8 +116,24 @@ class PETScVector
inline Vec& get() { return vec; }
protected:
std::string name_;
MPI_Comm petscCommunicator_;
Vec vec;
bool allocated_;
private:
inline void allocateVector( uint_t localSize )
{