Skip to content

GitLab

Commit ba5d0a3e authored by wagnandr's avatar wagnandr
Browse files

updates gmres apps

parent cd5bacdb
Hide whitespace changes
Inline Side-by-side
apps/StokesGMRES/StokesGMRES.cpp
View file @ ba5d0a3e
......@@ -25,7 +25,7 @@
#include "core/mpi/MPIManager.h"
#include "hyteg/composites/P1StokesFunction.hpp"
#include "hyteg/composites/P1StokesOperator.hpp"
#include "hyteg/composites/P1P1StokesOperator.hpp"
#include "hyteg/dataexport/VTKOutput.hpp"
#include "hyteg/functions/FunctionProperties.hpp"
#include "hyteg/gridtransferoperators/P1P1StokesToP1P1StokesProlongation.hpp"
......@@ -127,13 +127,13 @@ int main( int argc, char* argv[] )
hyteg::VTKOutput vtkOutput("./output", "StokesGMRES", storage);
if( mainConf.getParameter< bool >( "VTKOutput" ) )
{
vtkOutput.add( u.uvw );
vtkOutput.add( u.p );
vtkOutput.add( f.uvw );
vtkOutput.add( f.p );
vtkOutput.add( u.uvw() );
vtkOutput.add( u.p() );
vtkOutput.add( f.uvw() );
vtkOutput.add( f.p() );
}
hyteg::P1StokesOperator L( storage, minLevel, maxLevel );
hyteg::P1P1StokesOperator L( storage, minLevel, maxLevel );
std::function< real_t( const hyteg::Point3D& ) > rhsPlumeX = [sourcePoint, sourceRadius]( const hyteg::Point3D& x ) {
const real_t distToSourcePoint = ( x - sourcePoint ).norm();
......@@ -162,7 +162,7 @@ int main( int argc, char* argv[] )
std::function< real_t( const hyteg::Point3D& ) > zero = []( const hyteg::Point3D& ) { return 0.0; };
std::function< real_t( const hyteg::Point3D& ) > ones = []( const hyteg::Point3D& ) { return 1.0; };
f.uvw.interpolate( {rhsPlumeX, rhsPlumeY, rhsPlumeZ}, maxLevel );
f.uvw().interpolate( {rhsPlumeX, rhsPlumeY, rhsPlumeZ}, maxLevel );
if( mainConf.getParameter< bool >( "VTKOutput" ) )
{
......@@ -173,7 +173,7 @@ int main( int argc, char* argv[] )
if( solverType == "plainGMRES" )
{
typedef hyteg::GMRESSolver< hyteg::P1StokesOperator >GMRESsolv;
typedef hyteg::GMRESSolver< hyteg::P1P1StokesOperator >GMRESsolv;
auto KleinerFeinerGMRES = GMRESsolv( storage, minLevel, maxLevel, maxKrylowDim, restartLength, arnoldiTOL, approxTOL, doubleOrthoTOL);
KleinerFeinerGMRES.solve( L, u, f, maxLevel );
......@@ -184,11 +184,11 @@ int main( int argc, char* argv[] )
} else if( solverType == "simplePrecGMRES" )
{
/// A block Preconditioner for GMRES /////
typedef StokesBlockDiagonalPreconditioner< hyteg::P1StokesOperator, hyteg::P1LumpedInvMassOperator >Preconditioner_T;
typedef StokesBlockDiagonalPreconditioner< hyteg::P1P1StokesOperator, hyteg::P1LumpedInvMassOperator >Preconditioner_T;
auto prec = std::make_shared< Preconditioner_T >( storage, minLevel, maxLevel, 5 );
typedef hyteg::GMRESSolver< hyteg::P1StokesOperator >GMRESsolv;
typedef hyteg::GMRESSolver< hyteg::P1P1StokesOperator >GMRESsolv;
auto KleinerFeinerGMRES = GMRESsolv( storage, minLevel, maxLevel, maxKrylowDim, restartLength, arnoldiTOL, approxTOL, doubleOrthoTOL, prec );
KleinerFeinerGMRES.solve( L, u, f, maxLevel );
......@@ -198,7 +198,7 @@ int main( int argc, char* argv[] )
WALBERLA_LOG_INFO_ON_ROOT("currentResidualL2 : " << currentResidualL2);
} else if( solverType == "plainMINRES" )
{
typedef hyteg::MinResSolver< hyteg::P1StokesOperator >MinResSolv;
typedef hyteg::MinResSolver< hyteg::P1P1StokesOperator >MinResSolv;
auto KleinerFeinerMinRes = MinResSolv( storage, minLevel, maxLevel, maxKrylowDim, approxTOL );
KleinerFeinerMinRes.solve( L, u, f, maxLevel );
......@@ -209,11 +209,11 @@ int main( int argc, char* argv[] )
} else if( solverType == "simplePrecMINRES" )
{
/// A block Preconditioner for MinRes /////
typedef StokesBlockDiagonalPreconditioner< hyteg::P1StokesOperator, hyteg::P1LumpedInvMassOperator >Preconditioner_T;
typedef StokesBlockDiagonalPreconditioner< hyteg::P1P1StokesOperator, hyteg::P1LumpedInvMassOperator >Preconditioner_T;
auto prec = std::make_shared< Preconditioner_T >( storage, minLevel, maxLevel, 5 );
typedef hyteg::MinResSolver< hyteg::P1StokesOperator >MinResSolv;
typedef hyteg::MinResSolver< hyteg::P1P1StokesOperator >MinResSolv;
auto KleinerFeinerMinRes = MinResSolv( storage, minLevel, maxLevel, maxKrylowDim, approxTOL, prec );
KleinerFeinerMinRes.solve( L, u, f, maxLevel );
......@@ -223,21 +223,21 @@ int main( int argc, char* argv[] )
WALBERLA_LOG_INFO_ON_ROOT("currentResidualL2 : " << currentResidualL2);
} else if( solverType == "multigrid")
{
typedef StokesPressureBlockPreconditioner< hyteg::P1StokesOperator, hyteg::P1LumpedInvMassOperator >PressurePreconditioner_T;
typedef StokesPressureBlockPreconditioner< hyteg::P1P1StokesOperator, hyteg::P1LumpedInvMassOperator >PressurePreconditioner_T;
auto pressurePrec = std::make_shared< PressurePreconditioner_T>( storage, minLevel, minLevel );
typedef hyteg::GMRESSolver< hyteg::P1StokesOperator >PressurePreconditionedGMRES_T;
typedef hyteg::GMRESSolver< hyteg::P1P1StokesOperator >PressurePreconditionedGMRES_T;
auto pressurePreconditionedGMRESSolver = std::make_shared< PressurePreconditionedGMRES_T >(
storage, minLevel, maxLevel, maxKrylowDim, restartLength, arnoldiTOL, approxTOL, doubleOrthoTOL, pressurePrec );
typedef GeometricMultigridSolver< hyteg::P1StokesOperator >GMRESMultigrid_T;
typedef GeometricMultigridSolver< hyteg::P1P1StokesOperator >GMRESMultigrid_T;
auto stokesRestriction = std::make_shared< hyteg::P1P1StokesToP1P1StokesRestriction>();
auto stokesProlongation = std::make_shared< hyteg::P1P1StokesToP1P1StokesProlongation >();
auto gaussSeidel = std::make_shared< hyteg::GaussSeidelSmoother< hyteg::P1StokesOperator::VelocityOperator_T > >();
auto multigridVelocityPreconditioner = std::make_shared< hyteg::StokesVelocityBlockBlockDiagonalPreconditioner< hyteg::P1StokesOperator > >( storage, gaussSeidel );
auto gaussSeidel = std::make_shared< hyteg::GaussSeidelSmoother< hyteg::P1P1StokesOperator::VelocityOperator_T > >();
auto multigridVelocityPreconditioner = std::make_shared< hyteg::StokesVelocityBlockBlockDiagonalPreconditioner< hyteg::P1P1StokesOperator > >( storage, gaussSeidel );
auto GMRESMultigridSmoother = std::make_shared< hyteg::UzawaSmoother< P1StokesOperator > >(storage, multigridVelocityPreconditioner,minLevel, maxLevel, 0.3);
auto GMRESMultigridSmoother = std::make_shared< hyteg::UzawaSmoother< P1P1StokesOperator > >(storage, multigridVelocityPreconditioner,minLevel, maxLevel, 0.3);
GMRESMultigrid_T GMRESsolv(
storage, GMRESMultigridSmoother, pressurePreconditionedGMRESSolver, stokesRestriction, stokesProlongation, minLevel, maxLevel, 2, 2, 2 );
......@@ -282,7 +282,7 @@ int main( int argc, char* argv[] )
uint_t globalSize = 0;
const uint_t localSize = numerator->enumerate(level, globalSize);
PETScManager petscManager( &argc, &argv );
PETScLUSolver< real_t, hyteg::P1StokesFunction, hyteg::P1StokesOperator > petScLUSolver( numerator, localSize, globalSize );
PETScLUSolver< real_t, hyteg::P1StokesFunction, hyteg::P1P1StokesOperator > petScLUSolver( numerator, localSize, globalSize );
f.u.assign( {1.0}, {&u.u}, level, DirichletBoundary );
f.v.assign( {1.0}, {&u.v}, level, DirichletBoundary );
f.w.assign( {1.0}, {&u.w}, level, DirichletBoundary );
......
src/hyteg/solvers/GMRESSolver.hpp
View file @ ba5d0a3e
......@@ -51,7 +51,6 @@ class GMRESSolver : public Solver< OperatorType >
real_t arnoldiTOL = 1e-16,
real_t approxTOL = 1e-16,
real_t doubleOrthoTOL = 0,
hyteg::DoFType flag = hyteg::Inner | hyteg::NeumannBoundary | hyteg::FreeslipBoundary,
std::shared_ptr< Solver< OperatorType > > preconditioner = std::make_shared< IdentityPreconditioner< OperatorType > >()
)
: storage_( storage )
......@@ -63,7 +62,7 @@ class GMRESSolver : public Solver< OperatorType >
, approxTOL_( approxTOL )
, doubleOrthoTOL_( doubleOrthoTOL )
, preconditioner_( preconditioner )
, flag_( flag )
, flag_( hyteg::Inner | hyteg::NeumannBoundary | hyteg::FreeslipBoundary )
, printInfo_( false )
, r0_( "r0", storage_, minLevel_, maxLevel_ )
, rPrec_( "rPrec", storage_, minLevel_, maxLevel_ )
......@@ -77,7 +76,7 @@ class GMRESSolver : public Solver< OperatorType >
// GMRES ACCORDING TO Y.SAAD
void solve( const OperatorType& A, const FunctionType& x, const FunctionType& b, const uint_t level ) override
{
{
timingTree_->start( "GMRES Solver" );
double approxERR = approxTOL_ + 1;
......@@ -85,8 +84,8 @@ class GMRESSolver : public Solver< OperatorType >
init(A, x, b, level, true);
for(uint_t j = 1; j < maxKrylowDim_; j++)
{
for(uint_t j = 1; j < maxKrylowDim_; j++)
{
if(j % restartLength_ == 0 )
{
generateFinalApproximation( x, level );
......@@ -98,22 +97,22 @@ class GMRESSolver : public Solver< OperatorType >
int currentIndex = j % restartLength_;
if(isOnFirstRestart)
{
if(isOnFirstRestart)
{
FunctionType w( "w", storage_, minLevel_, maxLevel_ );
w.copyBoundaryConditionFromFunction( x );
vecV_.push_back( w );
} else
{
vecV_[currentIndex].interpolate( 0.0 , level, flag_ );
}
}
A.apply( vecV_[currentIndex - 1], vecV_[currentIndex], level, flag_ );
wPrec_.interpolate( 0.0, level, hyteg::Inner );
wPrec_.copyBoundaryConditionFromFunction( x );
preconditioner_->solve( A, wPrec_, vecV_[currentIndex], level );
vecV_[currentIndex].assign( {1.0}, {wPrec_}, level, flag_ );
H_ = matrixResize( H_, currentIndex + 1, currentIndex );
for(uint_t i = 1; i <= currentIndex; i++)
for(uint_t i = 1; i <= currentIndex; i++)
{
double h = vecV_[currentIndex].dotGlobal( vecV_[i - 1], level, flag_ );
H_(i-1, currentIndex - 1) = h;
......@@ -121,17 +120,17 @@ class GMRESSolver : public Solver< OperatorType >
}
A.apply( vecV_[currentIndex - 1], orthoDiff_, level, flag_ );
orthoDiff_.add( {-1.0}, {vecV_[currentIndex]}, level, flag_ );
if( std::sqrt(orthoDiff_.dotGlobal( orthoDiff_, level, flag_ )) > doubleOrthoTOL_ )
if( std::sqrt(orthoDiff_.dotGlobal( orthoDiff_, level, flag_ )) > doubleOrthoTOL_ )
{
if( printInfo_ ) {
WALBERLA_LOG_INFO_ON_ROOT(" [GMRES] invoked double-orthogonalization at iteration " << j );
}
for(uint_t i = 1; i <= (currentIndex); i++)
for(uint_t i = 1; i <= (currentIndex); i++)
{
double h = vecV_[currentIndex].dotGlobal( vecV_[i-1], level, flag_ );
H_(i-1, currentIndex - 1) += h;
vecV_[currentIndex].add( {-h}, {vecV_[i % restartLength_ - 1]}, level, flag_ );
}
}
}
double wNorm = std::sqrt( vecV_[currentIndex].dotGlobal( vecV_[currentIndex], level, flag_ ));
......@@ -142,7 +141,7 @@ class GMRESSolver : public Solver< OperatorType >
}
vecV_[currentIndex].assign( {1.0/wNorm}, {vecV_[currentIndex]}, level, flag_ );
if(wNorm <= arnoldiTOL_ || approxERR <= approxTOL_)
if(wNorm <= arnoldiTOL_ || approxERR <= approxTOL_)
{
break;
}
......@@ -154,7 +153,7 @@ class GMRESSolver : public Solver< OperatorType >
}
private:
void init(const OperatorType& A, const FunctionType& x, const FunctionType& b, const uint_t level, bool isFirstInit)
void init(const OperatorType& A, const FunctionType& x, const FunctionType& b, const uint_t level, bool isFirstInit)
{
A.apply(x, r0_, level, flag_);
r0_.assign( {1.0, -1.0}, {b, r0_}, level, flag_ );
......@@ -167,12 +166,12 @@ private:
wPrec_.copyBoundaryConditionFromFunction( x );
orthoDiff_.copyBoundaryConditionFromFunction( x );
beta_ = std::sqrt( r0_.dotGlobal( r0_, level, flag_ ) );
beta_ = std::sqrt( r0_.dotGlobal( r0_, level, flag_ ) );
H_ = Eigen::MatrixXd::Zero(0,0);
Q_ = Eigen::MatrixXd::Ones(1,1);
if(isFirstInit)
if(isFirstInit)
{
FunctionType v0( "v0", storage_, minLevel_, maxLevel_ );
v0.assign( {1.0/beta_}, {r0_}, level, flag_ );
......@@ -181,7 +180,7 @@ private:
{
vecV_[0].assign( {1.0/beta_}, {r0_}, level, flag_ );
}
}
}
const std::shared_ptr< PrimitiveStorage >& storage_;
uint_t minLevel_;
......@@ -207,22 +206,22 @@ private:
std::shared_ptr< walberla::WcTimingTree > timingTree_;
void generateFinalApproximation(const FunctionType& x, uint_t level)
void generateFinalApproximation(const FunctionType& x, uint_t level)
{
for(int i = 0; i < y_.size(); i++)
for(int i = 0; i < y_.size(); i++)
{
x.add( {y_(i)}, {vecV_[i]}, level, flag_ );
}
}
Eigen::VectorXd getUnitVector(int length, int j)
Eigen::VectorXd getUnitVector(int length, int j)
{
Eigen::VectorXd answer = Eigen::VectorXd::Zero(length);
answer(j) = 1;
return answer;
}
Eigen::VectorXd getUnitVector(int length, int j, double val)
Eigen::VectorXd getUnitVector(int length, int j, double val)
{
Eigen::VectorXd answer = Eigen::VectorXd::Zero(length);
answer(j) = val;
......@@ -231,7 +230,7 @@ private:
Eigen::MatrixXd getIdentity(int rowsNew, int colsNew) {
Eigen::MatrixXd answer = Eigen::MatrixXd::Zero(rowsNew, colsNew);
for(int i = 0; i < std::min(rowsNew, colsNew); i++)
for(int i = 0; i < std::min(rowsNew, colsNew); i++)
{
answer(i, i) = 1;
}
......@@ -241,9 +240,9 @@ private:
Eigen::MatrixXd matrixResize(Eigen::MatrixXd inputMatrix, int rowsNew, int colsNew)
{
Eigen::MatrixXd outputMatrix = Eigen::MatrixXd::Zero(rowsNew, colsNew);
for(uint_t j = 0; j < std::min(rowsNew, (int) inputMatrix.rows()); j++)
for(uint_t j = 0; j < std::min(rowsNew, (int) inputMatrix.rows()); j++)
{
for(int i = 0; i < std::min(colsNew, (int) inputMatrix.cols()); i++)
for(int i = 0; i < std::min(colsNew, (int) inputMatrix.cols()); i++)
{
outputMatrix(j, i) = inputMatrix(j, i);
}
......@@ -251,12 +250,12 @@ private:
return outputMatrix;
}
Eigen::MatrixXd expandQMatrix(Eigen::MatrixXd inputMatrix, int expandBy)
Eigen::MatrixXd expandQMatrix(Eigen::MatrixXd inputMatrix, int expandBy)
{
Eigen::MatrixXd outputMatrix = getIdentity(inputMatrix.rows() + expandBy, inputMatrix.cols() + expandBy);
for(int i = 0; i < inputMatrix.rows(); i++)
for(int i = 0; i < inputMatrix.rows(); i++)
{
for(int j = 0; j < inputMatrix.cols(); j++)
for(int j = 0; j < inputMatrix.cols(); j++)
{
outputMatrix(i,j) = inputMatrix(i,j);
}
......@@ -264,13 +263,13 @@ private:
return outputMatrix;
}
Eigen::VectorXd triangSolver(Eigen::MatrixXd triMat, Eigen::VectorXd targetVector)
Eigen::VectorXd triangSolver(Eigen::MatrixXd triMat, Eigen::VectorXd targetVector)
{
Eigen::VectorXd answer = Eigen::VectorXd::Zero(triMat.cols());
for(int i = triMat.cols()-1; i >= 0; i--)
for(int i = triMat.cols()-1; i >= 0; i--)
{
double targetOffset = 0;
for(int j = triMat.cols()-1; j > i; j--)
for(int j = triMat.cols()-1; j > i; j--)
{
targetOffset += triMat(i, j) * answer(j);
}
......@@ -279,23 +278,23 @@ private:
return answer;
}
Eigen::VectorXd hessenbergMinimizer(double beta, Eigen::MatrixXd& H, Eigen::MatrixXd& Q, int numUnfinishedColumns, double& approxERR)
Eigen::VectorXd hessenbergMinimizer(double beta, Eigen::MatrixXd& H, Eigen::MatrixXd& Q, int numUnfinishedColumns, double& approxERR)
{
Eigen::VectorXd approxVector;
if(H.rows() != Q.rows())
if(H.rows() != Q.rows())
{
Q = expandQMatrix(Q, H.rows() - Q.rows());
}
for(int i = H.cols() - numUnfinishedColumns; i < H.cols(); i++)
for(int i = H.cols() - numUnfinishedColumns; i < H.cols(); i++)
{
H.col(i) = Q * H.col(i);
}
for(int i = H.cols() - numUnfinishedColumns; i < H.cols(); i++)
for(int i = H.cols() - numUnfinishedColumns; i < H.cols(); i++)
{
if(H(i+1, i) == 0)
if(H(i+1, i) == 0)
{
continue;
}
......@@ -311,7 +310,7 @@ private:
Q = QNew * Q;
}
Eigen::MatrixXd equationLeftSide = H.block(0, 0, H.cols(), H.cols());
Eigen::MatrixXd equationLeftSide = H.block(0, 0, H.cols(), H.cols());
Eigen::VectorXd targetVector = Q * getUnitVector(H.rows(), 0, beta);
Eigen::VectorXd equationRightSide = targetVector.head(H.cols());
approxERR = std::abs(targetVector(targetVector.rows() - 1));
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment