Commit b0c4dcb2 authored by Andreas Burkhart's avatar Andreas Burkhart
Browse files

more changes to conv test

parent b54b504f
Pipeline #41173 passed with stages
in 169 minutes and 26 seconds
......@@ -19,187 +19,192 @@
*/
#pragma once
#include "core/timing/TimingTree.h"
#include "core/Format.hpp"
#include "core/timing/TimingTree.h"
#include "hyteg/solvers/preconditioners/IdentityPreconditioner.hpp"
#include "hyteg/solvers/Solver.hpp"
#include "hyteg/solvers/preconditioners/IdentityPreconditioner.hpp"
namespace hyteg {
template<class OperatorType >
template < class OperatorType >
class MinResSolver : public Solver< OperatorType >
{
public:
typedef typename OperatorType::srcType FunctionType;
MinResSolver(
const std::shared_ptr< PrimitiveStorage >& storage,
size_t minLevel,
size_t maxLevel,
uint_t maxIter = std::numeric_limits< uint_t >::max(),
real_t tolerance = 1e-16,
std::shared_ptr< Solver< OperatorType > > preconditioner = std::make_shared< IdentityPreconditioner< OperatorType > >() )
: maxIter_( maxIter )
, relativeTolerance_( tolerance )
, absoluteTolerance_( 1e-16 )
, printInfo_( false )
, flag_( hyteg::Inner | hyteg::NeumannBoundary | hyteg::FreeslipBoundary )
, preconditioner_( preconditioner )
, p_vm( "minres_vm", storage, minLevel, maxLevel )
, p_v( "minres_v", storage, minLevel, maxLevel )
, p_vp( "minres_vp", storage, minLevel, maxLevel )
, p_z( "minres_z", storage, minLevel, maxLevel )
, p_zp( "minres_zp", storage, minLevel, maxLevel )
, p_wm( "minres_wm", storage, minLevel, maxLevel )
, p_w( "minres_w", storage, minLevel, maxLevel )
, p_wp( "minres_wp", storage, minLevel, maxLevel )
, p_tmp( "minres_tmp", storage, minLevel, maxLevel )
, r_( "minres_r", storage, minLevel, maxLevel )
, timingTree_( storage->getTimingTree() )
{}
void solve( const OperatorType& A,const FunctionType& x, const FunctionType& b, const uint_t level ) override
{
timingTree_->start( "MinRes Solver" );
p_vm.copyBoundaryConditionFromFunction( x );
p_v.copyBoundaryConditionFromFunction( x );
p_vp.copyBoundaryConditionFromFunction( x );
p_z.copyBoundaryConditionFromFunction( x );
p_zp.copyBoundaryConditionFromFunction( x );
p_wm.copyBoundaryConditionFromFunction( x );
p_w.copyBoundaryConditionFromFunction( x );
p_wp.copyBoundaryConditionFromFunction( x );
p_tmp.copyBoundaryConditionFromFunction( x );
r_.copyBoundaryConditionFromFunction( x );
std::function<real_t(const hyteg::Point3D&)> zero = [](const hyteg::Point3D&) { return 0.0; };
p_vm.interpolate(zero, level, flag_);
p_wm.interpolate(zero, level, flag_);
p_w.interpolate(zero, level, flag_);
A.apply(x, r_, level, flag_);
p_v.assign({1.0, -1.0}, {b, r_}, level, flag_);
p_z.interpolate(0, level, flag_);
preconditioner_->solve(A, p_z, p_v, level );
real_t gamma_old = 1.0;
real_t gamma_new = std::sqrt(p_z.dotGlobal(p_v, level, flag_));
real_t res_start = gamma_new;
real_t eta = gamma_new;
real_t s_old = 0.0;
real_t s_new = 0.0;
real_t c_old = 1.0;
real_t c_new = 1.0;
if (printInfo_) {
WALBERLA_LOG_INFO_ON_ROOT("[MinRes] residuum: "<< std::scientific << std::abs(gamma_new));
}
if (gamma_new < absoluteTolerance_ )
{
if (printInfo_) {
WALBERLA_LOG_INFO_ON_ROOT("[MinRes] converged");
public:
typedef typename OperatorType::srcType FunctionType;
MinResSolver(
const std::shared_ptr< PrimitiveStorage >& storage,
size_t minLevel,
size_t maxLevel,
uint_t maxIter = std::numeric_limits< uint_t >::max(),
real_t tolerance = 1e-16,
std::shared_ptr< Solver< OperatorType > > preconditioner = std::make_shared< IdentityPreconditioner< OperatorType > >() )
: maxIter_( maxIter )
, relativeTolerance_( tolerance )
, absoluteTolerance_( 1e-16 )
, printInfo_( false )
, flag_( hyteg::Inner | hyteg::NeumannBoundary | hyteg::FreeslipBoundary )
, preconditioner_( preconditioner )
, p_vm( "minres_vm", storage, minLevel, maxLevel )
, p_v( "minres_v", storage, minLevel, maxLevel )
, p_vp( "minres_vp", storage, minLevel, maxLevel )
, p_z( "minres_z", storage, minLevel, maxLevel )
, p_zp( "minres_zp", storage, minLevel, maxLevel )
, p_wm( "minres_wm", storage, minLevel, maxLevel )
, p_w( "minres_w", storage, minLevel, maxLevel )
, p_wp( "minres_wp", storage, minLevel, maxLevel )
, p_tmp( "minres_tmp", storage, minLevel, maxLevel )
, r_( "minres_r", storage, minLevel, maxLevel )
, timingTree_( storage->getTimingTree() )
, numberOfIterations_( maxIter )
{}
void solve( const OperatorType& A, const FunctionType& x, const FunctionType& b, const uint_t level ) override
{
timingTree_->start( "MinRes Solver" );
p_vm.copyBoundaryConditionFromFunction( x );
p_v.copyBoundaryConditionFromFunction( x );
p_vp.copyBoundaryConditionFromFunction( x );
p_z.copyBoundaryConditionFromFunction( x );
p_zp.copyBoundaryConditionFromFunction( x );
p_wm.copyBoundaryConditionFromFunction( x );
p_w.copyBoundaryConditionFromFunction( x );
p_wp.copyBoundaryConditionFromFunction( x );
p_tmp.copyBoundaryConditionFromFunction( x );
r_.copyBoundaryConditionFromFunction( x );
std::function< real_t( const hyteg::Point3D& ) > zero = []( const hyteg::Point3D& ) { return 0.0; };
p_vm.interpolate( zero, level, flag_ );
p_wm.interpolate( zero, level, flag_ );
p_w.interpolate( zero, level, flag_ );
A.apply( x, r_, level, flag_ );
p_v.assign( { 1.0, -1.0 }, { b, r_ }, level, flag_ );
p_z.interpolate( 0, level, flag_ );
preconditioner_->solve( A, p_z, p_v, level );
real_t gamma_old = 1.0;
real_t gamma_new = std::sqrt( p_z.dotGlobal( p_v, level, flag_ ) );
real_t res_start = gamma_new;
real_t eta = gamma_new;
real_t s_old = 0.0;
real_t s_new = 0.0;
real_t c_old = 1.0;
real_t c_new = 1.0;
if ( printInfo_ )
{
WALBERLA_LOG_INFO_ON_ROOT( "[MinRes] residuum: " << std::scientific << std::abs( gamma_new ) );
}
timingTree_->stop( "MinRes Solver" );
return;
}
for(size_t i = 0; i < maxIter_; ++i) {
p_z.assign({real_t(1) / gamma_new}, {p_z}, level, flag_);
A.apply(p_z, p_vp, level, flag_);
real_t delta = p_vp.dotGlobal(p_z, level, flag_);
p_vp.assign({1.0, -delta / gamma_new, -gamma_new / gamma_old}, {p_vp, p_v, p_vm}, level, flag_);
p_zp.interpolate(0, level, flag_);
preconditioner_->solve(A, p_zp, p_vp, level );
gamma_old = gamma_new;
gamma_new = std::sqrt(p_zp.dotGlobal(p_vp, level, flag_));
real_t alpha0 = c_new * delta - c_old * s_new * gamma_old;
real_t alpha1 = std::sqrt(alpha0 * alpha0 + gamma_new * gamma_new);
real_t alpha2 = s_new * delta + c_old * c_new * gamma_old;
real_t alpha3 = s_old * gamma_old;
c_old = c_new;
c_new = alpha0 / alpha1;
s_old = s_new;
s_new = gamma_new / alpha1;
p_wp.assign({real_t(1)/alpha1, -alpha3/alpha1, -alpha2/alpha1}, {p_z, p_wm, p_w}, level, flag_);
x.add({c_new * eta}, {p_wp}, level, flag_);
eta = -s_new * eta;
p_tmp.swap( p_vp, level );
p_vp.swap( p_vm, level );
p_vm.swap( p_v, level );
p_v.swap( p_tmp, level );
p_tmp.swap( p_wp, level );
p_wp.swap( p_wm, level );
p_wm.swap( p_w, level );
p_w.swap( p_tmp, level );
p_tmp.swap( p_zp, level );
p_zp.swap( p_z, level );
p_z.swap( p_tmp, level );
if (printInfo_)
if ( gamma_new < absoluteTolerance_ )
{
WALBERLA_LOG_INFO_ON_ROOT(walberla::format("[MinRes] iter: %6d | residuum: %10.5e", i, std::abs(eta)));
if ( printInfo_ )
{
WALBERLA_LOG_INFO_ON_ROOT( "[MinRes] converged" );
}
timingTree_->stop( "MinRes Solver" );
return;
}
if (std::abs(eta)/res_start < relativeTolerance_ || std::abs(eta) < absoluteTolerance_ )
for ( size_t i = 0; i < maxIter_; ++i )
{
if (printInfo_)
{
WALBERLA_LOG_INFO_ON_ROOT("[MinRes] converged after " << std::defaultfloat << i << " iterations");
}
break;
p_z.assign( { real_t( 1 ) / gamma_new }, { p_z }, level, flag_ );
A.apply( p_z, p_vp, level, flag_ );
real_t delta = p_vp.dotGlobal( p_z, level, flag_ );
p_vp.assign( { 1.0, -delta / gamma_new, -gamma_new / gamma_old }, { p_vp, p_v, p_vm }, level, flag_ );
p_zp.interpolate( 0, level, flag_ );
preconditioner_->solve( A, p_zp, p_vp, level );
gamma_old = gamma_new;
gamma_new = std::sqrt( p_zp.dotGlobal( p_vp, level, flag_ ) );
real_t alpha0 = c_new * delta - c_old * s_new * gamma_old;
real_t alpha1 = std::sqrt( alpha0 * alpha0 + gamma_new * gamma_new );
real_t alpha2 = s_new * delta + c_old * c_new * gamma_old;
real_t alpha3 = s_old * gamma_old;
c_old = c_new;
c_new = alpha0 / alpha1;
s_old = s_new;
s_new = gamma_new / alpha1;
p_wp.assign( { real_t( 1 ) / alpha1, -alpha3 / alpha1, -alpha2 / alpha1 }, { p_z, p_wm, p_w }, level, flag_ );
x.add( { c_new * eta }, { p_wp }, level, flag_ );
eta = -s_new * eta;
p_tmp.swap( p_vp, level );
p_vp.swap( p_vm, level );
p_vm.swap( p_v, level );
p_v.swap( p_tmp, level );
p_tmp.swap( p_wp, level );
p_wp.swap( p_wm, level );
p_wm.swap( p_w, level );
p_w.swap( p_tmp, level );
p_tmp.swap( p_zp, level );
p_zp.swap( p_z, level );
p_z.swap( p_tmp, level );
if ( printInfo_ )
{
WALBERLA_LOG_INFO_ON_ROOT( walberla::format( "[MinRes] iter: %6d | residuum: %10.5e", i, std::abs( eta ) ) );
}
if ( std::abs( eta ) / res_start < relativeTolerance_ || std::abs( eta ) < absoluteTolerance_ )
{
numberOfIterations_ = i;
if ( printInfo_ )
{
WALBERLA_LOG_INFO_ON_ROOT( "[MinRes] converged after " << std::defaultfloat << i << " iterations" );
}
break;
}
}
}
timingTree_->stop( "MinRes Solver" );
}
void setPrintInfo( bool printInfo ) { printInfo_ = printInfo; }
void setAbsoluteTolerance( real_t absoluteTolerance ) { absoluteTolerance_ = absoluteTolerance; }
void setRelativeTolerance( real_t relativeTolerance ) { relativeTolerance_ = relativeTolerance; }
timingTree_->stop( "MinRes Solver" );
}
private:
uint_t getIterations() { return numberOfIterations_; }
void setPrintInfo( bool printInfo ) { printInfo_ = printInfo; }
void setAbsoluteTolerance( real_t absoluteTolerance ) { absoluteTolerance_ = absoluteTolerance; }
void setRelativeTolerance( real_t relativeTolerance ) { relativeTolerance_ = relativeTolerance; }
uint_t maxIter_;
real_t relativeTolerance_;
real_t absoluteTolerance_;
bool printInfo_;
hyteg::DoFType flag_;
private:
uint_t maxIter_;
real_t relativeTolerance_;
real_t absoluteTolerance_;
bool printInfo_;
hyteg::DoFType flag_;
uint_t numberOfIterations_;
std::shared_ptr< Solver< OperatorType > > preconditioner_;
std::shared_ptr< Solver< OperatorType > > preconditioner_;
FunctionType p_vm;
FunctionType p_v;
FunctionType p_vp;
FunctionType p_vm;
FunctionType p_v;
FunctionType p_vp;
FunctionType p_z;
FunctionType p_zp;
FunctionType p_z;
FunctionType p_zp;
FunctionType p_wm;
FunctionType p_w;
FunctionType p_wp;
FunctionType p_wm;
FunctionType p_w;
FunctionType p_wp;
FunctionType p_tmp;
FunctionType p_tmp;
FunctionType r_;
FunctionType r_;
std::shared_ptr< walberla::WcTimingTree > timingTree_;
std::shared_ptr< walberla::WcTimingTree > timingTree_;
};
}
} // namespace hyteg
......@@ -50,6 +50,7 @@
#include "hyteg/solvers/GMRESSolver.hpp"
#include "hyteg/solvers/GaussSeidelSmoother.hpp"
#include "hyteg/solvers/MinresSolver.hpp"
#include "hyteg/solvers/WeightedJacobiSmoother.hpp"
#include "hyteg/solvers/preconditioners/stokes/StokesBlockDiagonalPreconditioner.hpp"
#include "hyteg/solvers/preconditioners/stokes/StokesPressureBlockPreconditioner.hpp"
#include "hyteg/solvers/preconditioners/stokes/StokesVelocityBlockBlockDiagonalPreconditioner.hpp"
......@@ -64,26 +65,28 @@ using StokesOp = hyteg::P2P1TaylorHoodCompStokesOperator;
// parameters
const uint_t minLevel = 1;
const uint_t maxLevel = 4;
const uint_t maxLevel = 6;
const uint_t startingLevel = 1;
//const uint_t max_outer_iter = 5;
const uint_t max_coarse_iter = 200;
const real_t coarse_tolerance = 1e-16;
const real_t coarse_tolerance = 1e-12;
const real_t boundary_tolerance = 1e-16;
const uint_t maxKrylowDim = 500;
const uint_t restartLength = 200;
const real_t arnoldiTOL = 1e-16;
const real_t approxTOL = 1e-14;
const real_t approxTOL = 1e-10;
const real_t doubleOrthoTOL = 0;
const bool printInfo = true;
const bool printInfo = false;
const bool vtkOut = false;
// 0 = GMRES, 1 = MINRES
void check( const std::array< std::function< real_t( const hyteg::Point3D& ) >, 11 >& functions,
const hyteg::MeshInfo& meshInfo,
uint_t dim,
uint_t scenarioNumber = 0 )
uint_t scenarioNumber = 0,
uint_t solverType = 0 )
{
SetupPrimitiveStorage setupStorage( meshInfo, uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) );
......@@ -132,13 +135,16 @@ void check( const std::array< std::function< real_t( const hyteg::Point3D& ) >,
auto laplRestriction = std::make_shared< hyteg::P2toP2QuadraticRestriction >();
auto laplProlongation = std::make_shared< hyteg::P2toP2QuadraticProlongation >();
auto gaussSeidel = std::make_shared< hyteg::GaussSeidelSmoother< StokesOp::VelocityOperator_T > >();
auto coarseGridSolver = std::make_shared< hyteg::CGSolver< StokesOp::VelocityOperator_T > >(
storage, minLevel, minLevel, max_coarse_iter, coarse_tolerance );
//coarseGridSolver->setPrintInfo( true );
auto multigridVelocityPreconditioner = std::make_shared< hyteg::GeometricMultigridSolver< StokesOp::VelocityOperator_T > >(
storage, gaussSeidel, coarseGridSolver, laplRestriction, laplProlongation, minLevel, maxLevel );
storage, gaussSeidel, coarseGridSolver, laplRestriction, laplProlongation, minLevel, maxLevel, 5, 5 );
auto prec = std::make_shared< Preconditioner_T >( storage, minLevel, maxLevel, 3, multigridVelocityPreconditioner );
auto prec = std::make_shared< Preconditioner_T >( storage, minLevel, maxLevel, 1, multigridVelocityPreconditioner );
// P2P1 Prolong
auto stokesProlongation = std::make_shared< hyteg::P2P1StokesToP2P1StokesProlongation >();
......@@ -195,11 +201,25 @@ void check( const std::array< std::function< real_t( const hyteg::Point3D& ) >,
uint_t globalDoFsWorking = hyteg::numberOfGlobalDoFs< hyteg::P2P1TaylorHoodFunctionTag >( *storage, workingLevel );
WALBERLA_LOG_INFO_ON_ROOT( "--- Level " << workingLevel << " " << globalDoFsWorking << " DoFs ---" );
// define solver
auto solver = hyteg::GMRESSolver< StokesOp >(
storage, minLevel, maxLevel, maxKrylowDim, restartLength, arnoldiTOL, approxTOL, doubleOrthoTOL, prec );
std::shared_ptr< hyteg::Solver< StokesOp > > solver;
if ( solverType == 0 )
{
auto newsolver = std::make_shared< hyteg::GMRESSolver< StokesOp > >(
storage, minLevel, maxLevel, maxKrylowDim, restartLength, arnoldiTOL, approxTOL, doubleOrthoTOL, prec );
newsolver->setPrintInfo( printInfo );
solver = std::shared_ptr< hyteg::Solver< StokesOp > >( newsolver );
}
else
{
auto newsolver =
std::make_shared< MinResSolver< StokesOp > >( storage, minLevel, maxLevel, maxKrylowDim, approxTOL, prec );
newsolver->setPrintInfo( printInfo );
solver = std::shared_ptr< hyteg::Solver< StokesOp > >( newsolver );
}
//auto solver = MinResSolver< StokesOp >( storage, minLevel, maxLevel, maxKrylowDim );
solver.setPrintInfo( printInfo );
//solver->setPrintInfo( printInfo );
// real_t disc_l2_up = std::sqrt( up.dotGlobal( up, workingLevel ) / std::sqrt( (real_t) globalDoFsWorking ) );
// WALBERLA_LOG_INFO_ON_ROOT( "Disc L2 uvw: " << disc_l2_up );
......@@ -220,10 +240,24 @@ void check( const std::array< std::function< real_t( const hyteg::Point3D& ) >,
massP2.apply( f.uvw(), f2.uvw(), workingLevel, hyteg::Inner | hyteg::NeumannBoundary );
massP1.apply( f.p(), f2.p(), workingLevel, hyteg::Inner | hyteg::NeumannBoundary );
solver.solve( Op, up, f2, workingLevel );
solver->solve( Op, up, f2, workingLevel );
WALBERLA_LOG_INFO_ON_ROOT( "Number of GMRES iterations " << solver.getIterations() << "/" << maxKrylowDim
<< ", restarted every " << restartLength << " iterations" );
uint_t niter = 0;
if ( solverType == 0 )
{
std::shared_ptr< hyteg::GMRESSolver< StokesOp > > newsolver =
std::dynamic_pointer_cast< hyteg::GMRESSolver< StokesOp > >( solver );
niter = newsolver->getIterations();
}
else
{
std::shared_ptr< hyteg::MinResSolver< StokesOp > > newsolver =
std::dynamic_pointer_cast< hyteg::MinResSolver< StokesOp > >( solver );
niter = newsolver->getIterations();
}
WALBERLA_LOG_INFO_ON_ROOT( "Number of iterations " << niter << "/" << maxKrylowDim << ", restarted every " << restartLength
<< " iterations" );
hyteg::vertexdof::projectMean( up.p(), workingLevel );
......@@ -264,6 +298,9 @@ void check( const std::array< std::function< real_t( const hyteg::Point3D& ) >,
WALBERLA_LOG_INFO_ON_ROOT( "ratio_p = " << ratio_p );
}
//auto tree = storage->getTimingTree();
//WALBERLA_LOG_INFO_ON_ROOT( *tree );
if ( vtkOut )
{
vtkOutput.add( up.uvw()[0] );
......
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