Commit 7d66ddc1 authored by Marcus Mohr's avatar Marcus Mohr
Browse files

Refactors StokesPCGSolverOld to become dimension independent

Commit changes StokesPCGSolverOld so that it will work in 2D and 3D.
By using the new VectorFunctions and Div/DivT operators the solver
becomes agnostic implementationwise w.r.t. dimension.

Commit also enhances the P2P1StokesSchurCGConvergenceTest by adding a
test case for 3D. Convergence not fast so we currently just make five
iterations to see that it works technically.

Also renames P2P1StokesSchurCGConvergenceTest. New name is closer to
the class name: P2P1StokesPCGSolverConvergenceTest
parent e72d5a47
Pipeline #27287 passed with stages
in 95 minutes and 54 seconds
......@@ -80,8 +80,12 @@ class StokesPCGSolverOld : public Solver< OperatorType >
invMass( storage, minLevel, maxLevel )
{
const bool available = !storage->hasGlobalCells() && std::is_same< OperatorType, P2P1TaylorHoodStokesOperator >::value;
WALBERLA_CHECK( available, "SchurCG solver currently only available for Taylor-Hood discretization on 2D domains." );
const bool opCheck = std::is_same< OperatorType, P2P1TaylorHoodStokesOperator >::value;
WALBERLA_CHECK( opCheck, "StokesPCG solver currently only available for P2-P1-Taylor-Hood discretization." );
if ( storage->hasGlobalCells() )
{
WALBERLA_LOG_WARNING_ON_ROOT( "StokesPCG solver not fully tested for 3D settings!" );
}
}
void solve( const OperatorType& stokesOperator, const FunctionType& u, const FunctionType& b, const uint_t level )
......@@ -94,19 +98,15 @@ class StokesPCGSolverOld : public Solver< OperatorType >
/////////
// velocity
R.uvw.u.assign( {1.0}, {b.uvw.u}, level, flag_ );
R.uvw.v.assign( {1.0}, {b.uvw.v}, level, flag_ );
stokesOperator.A.apply( u.uvw.u, tmp.uvw.u, level, flag_ );
stokesOperator.A.apply( u.uvw.v, tmp.uvw.v, level, flag_ );
// stokesOperator.divT_x.apply( u.p, tmp.uvw.u, level, flag_, Add );
// stokesOperator.divT_y.apply( u.p, tmp.uvw.v, level, flag_, Add );
R.uvw.assign( {1.0}, {b.uvw}, level, flag_ );
for ( uint_t k = 0; k < u.uvw.getDimension(); k++ )
{
stokesOperator.A.apply( u.uvw[k], tmp.uvw[k], level, flag_ );
}
stokesOperator.divT.apply( u.p, tmp.uvw, level, flag_, Add );
R.uvw.u.add( {-1.0}, {tmp.uvw.u}, level, flag_ );
R.uvw.v.add( {-1.0}, {tmp.uvw.v}, level, flag_ );
R.uvw.add( {-1.0}, {tmp.uvw}, level, flag_ );
// pressure
// stokesOperator.div_x.apply( u.uvw.u, tmp.p, level, flag_ | DirichletBoundary, Replace );
// stokesOperator.div_y.apply( u.uvw.v, tmp.p, level, flag_ | DirichletBoundary, Add );
stokesOperator.div.apply( u.uvw, tmp.p, level, flag_ | DirichletBoundary, Replace );
R.p.assign( {-1.0}, {tmp.p}, level, flag_ | DirichletBoundary );
......@@ -115,8 +115,10 @@ class StokesPCGSolverOld : public Solver< OperatorType >
/////////////
// velocity
velocityBlockSolver_->solve( stokesOperator.A, RHat.uvw.u, R.uvw.u, level );
velocityBlockSolver_->solve( stokesOperator.A, RHat.uvw.v, R.uvw.v, level );
for ( uint_t k = 0; k < u.uvw.getDimension(); k++ )
{
velocityBlockSolver_->solve( stokesOperator.A, RHat.uvw[k], R.uvw[k], level );
}
// pressure
stokesOperator.div.apply( RHat.uvw, RHat.p, level, flag_ | DirichletBoundary, Replace );
......@@ -142,15 +144,17 @@ class StokesPCGSolverOld : public Solver< OperatorType >
alpha_n = RTilde.p.dotGlobal( RHat.p, level, flag_ | DirichletBoundary );
// calculate M * P
// 1. calc tmp := A^-1 * B^T * P_pressure
// stokesOperator.divT_x.apply( P.p, tmp.uvw.u, level, flag_ );
// stokesOperator.divT_y.apply( P.p, tmp.uvw.v, level, flag_ );
stokesOperator.divT.apply( P.p, tmp.uvw, level, flag_ );
velocityBlockSolver_->solve( stokesOperator.A, tmp_2.uvw.u, tmp.uvw.u, level );
velocityBlockSolver_->solve( stokesOperator.A, tmp_2.uvw.v, tmp.uvw.v, level );
for ( uint_t k = 0; k < tmp_2.uvw.getDimension(); k++ )
{
velocityBlockSolver_->solve( stokesOperator.A, tmp_2.uvw[k], tmp.uvw[k], level );
}
// 2. velocity of M * P
MP.uvw.u.assign( {1.0, 1.0}, {P.uvw.u, tmp_2.uvw.u}, level, flag_ );
MP.uvw.v.assign( {1.0, 1.0}, {P.uvw.v, tmp_2.uvw.v}, level, flag_ );
MP.uvw.assign( {1.0, 1.0}, {P.uvw, tmp_2.uvw}, level, flag_ );
// 3. pressure of M * P
stokesOperator.div.apply( tmp_2.uvw, MP.p, level, flag_ | DirichletBoundary, Replace );
......@@ -164,9 +168,7 @@ class StokesPCGSolverOld : public Solver< OperatorType >
// X //
///////
// u.add( {alpha}, {P}, level, flag_ );
u.uvw.u.add( {alpha}, {P.uvw.u}, level, flag_ );
u.uvw.v.add( {alpha}, {P.uvw.v}, level, flag_ );
u.uvw.add( {alpha}, {P.uvw}, level, flag_ );
u.p.add( {alpha}, {P.p}, level, flag_ | DirichletBoundary );
//////////////
......@@ -198,7 +200,7 @@ class StokesPCGSolverOld : public Solver< OperatorType >
stokesOperator.apply( u, Au, level, flag_ );
residual.assign( {1.0, -1.0}, {b, Au}, level, flag_ );
real_t currentResidual = std::sqrt( residual.dotGlobal( residual, level, flag_ ) ) / numGlobalDoFs;
WALBERLA_LOG_DEVEL_ON_ROOT( "CURRENT SCHUR CG RESIDUAL (" << i << "): " << std::scientific << currentResidual );
WALBERLA_LOG_DEVEL_ON_ROOT( "CURRENT PCG RESIDUAL (" << i << "): " << std::scientific << currentResidual );
if ( currentResidual < residualTolerance_ )
return;
......@@ -210,12 +212,13 @@ class StokesPCGSolverOld : public Solver< OperatorType >
if ( !solveExactly_ )
{
stokesOperator.A.apply( RHat.uvw.u, tmp.uvw.u, level, flag_ );
stokesOperator.A.apply( RHat.uvw.v, tmp.uvw.v, level, flag_ );
tmp.uvw.u.add( {-1.0}, {R.uvw.u}, level );
tmp.uvw.v.add( {-1.0}, {R.uvw.v}, level );
beta_n += RHat.uvw.u.dotGlobal( tmp.uvw.u, level, flag_ );
beta_n += RHat.uvw.v.dotGlobal( tmp.uvw.v, level, flag_ );
for ( uint_t k = 0; k < u.uvw.getDimension(); k++ )
{
stokesOperator.A.apply( RHat.uvw[k], tmp.uvw[k], level, flag_ );
}
tmp.uvw.add( {-1.0}, {R.uvw}, level );
beta_n += RHat.uvw.dotGlobal( tmp.uvw, level, flag_ );
}
beta_d = alpha_n;
......@@ -251,18 +254,18 @@ class StokesPCGSolverOld : public Solver< OperatorType >
alpha_n = beta_n;
// calculate M * P
// 1. calc tmp := A^-1 * B^T * P_pressure
// stokesOperator.divT_x.apply( P.p, tmp.uvw.u, level, flag_ );
// stokesOperator.divT_y.apply( P.p, tmp.uvw.v, level, flag_ );
stokesOperator.divT.apply( P.p, tmp.uvw, level, flag_ );
velocityBlockSolver_->solve( stokesOperator.A, tmp_2.uvw.u, tmp.uvw.u, level );
velocityBlockSolver_->solve( stokesOperator.A, tmp_2.uvw.v, tmp.uvw.v, level );
for ( uint_t k = 0; k < u.uvw.getDimension(); k++ )
{
velocityBlockSolver_->solve( stokesOperator.A, tmp_2.uvw[k], tmp.uvw[k], level );
}
// 2. velocity of M * P
MP.uvw.u.assign( {1.0, 1.0}, {P.uvw.u, tmp_2.uvw.u}, level, flag_ );
MP.uvw.v.assign( {1.0, 1.0}, {P.uvw.v, tmp_2.uvw.v}, level, flag_ );
MP.uvw.assign( {1.0, 1.0}, {P.uvw, tmp_2.uvw}, level, flag_ );
// 3. pressure of M * P
// stokesOperator.div_x.apply( tmp_2.uvw.u, MP.p, level, flag_ | DirichletBoundary, Replace );
// stokesOperator.div_y.apply( tmp_2.uvw.v, MP.p, level, flag_ | DirichletBoundary, Add );
stokesOperator.div.apply( tmp_2.uvw, MP.p, level, flag_ | DirichletBoundary, Replace );
alpha_d = P.p.dotGlobal( MP.p, level, flag_ | DirichletBoundary );
......@@ -291,8 +294,7 @@ class StokesPCGSolverOld : public Solver< OperatorType >
// X //
///////
u.uvw.u.add( {alpha}, {P.uvw.u}, level, flag_ );
u.uvw.v.add( {alpha}, {P.uvw.v}, level, flag_ );
u.uvw.add( {alpha}, {P.uvw}, level, flag_ );
u.p.add( {alpha}, {P.p}, level, flag_ | DirichletBoundary );
//////////////
......
......@@ -113,8 +113,8 @@ waLBerla_compile_test(FILES convergence/P2P1StokesMinResConvergenceTest.cpp DEPE
waLBerla_execute_test(NAME P2P1StokesMinResConvergenceTest)
if(CMAKE_BUILD_TYPE MATCHES "Release")
waLBerla_compile_test(FILES convergence/P2P1StokesSchurCGConvergenceTest.cpp DEPENDS hyteg core)
waLBerla_execute_test(NAME P2P1StokesSchurCGConvergenceTest)
waLBerla_compile_test(FILES convergence/P2P1StokesPCGSolverConvergenceTest.cpp DEPENDS hyteg core)
waLBerla_execute_test(NAME P2P1StokesPCGSolverConvergenceTest)
endif()
waLBerla_compile_test(FILES convergence/P2GSConvergenceTest.cpp DEPENDS hyteg core )
......
......@@ -41,9 +41,12 @@ using walberla::real_t;
using walberla::uint_c;
using walberla::uint_t;
using std::sin;
using std::cos;
namespace hyteg {
void P2P1SchurCGConvergenceTest( const uint_t & level, const MeshInfo & meshInfo )
void P2P1PCGSolverConvergenceTest( const uint_t & level, const MeshInfo & meshInfo )
{
SetupPrimitiveStorage setupStorage( meshInfo, uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) );
......@@ -54,7 +57,7 @@ void P2P1SchurCGConvergenceTest( const uint_t & level, const MeshInfo & meshInfo
hyteg::loadbalancing::roundRobin( setupStorage );
std::shared_ptr< PrimitiveStorage > storage = std::make_shared< PrimitiveStorage >( setupStorage );
writeDomainPartitioningVTK( storage, "../../output", "P2P1Stokes2DSchurCGConvergence_Domain" );
// writeDomainPartitioningVTK( storage, "../../output", "P2P1Stokes2DPCGSolverConvergence_Domain" );
hyteg::P2P1TaylorHoodFunction< real_t > x( "x", storage, minLevel, level );
hyteg::P2P1TaylorHoodFunction< real_t > x_exact( "x_exact", storage, minLevel, level );
......@@ -65,32 +68,58 @@ void P2P1SchurCGConvergenceTest( const uint_t & level, const MeshInfo & meshInfo
hyteg::P2P1TaylorHoodStokesOperator A( storage, minLevel, level );
std::function< real_t( const hyteg::Point3D& ) > exactU = []( const hyteg::Point3D& xx ) { return real_c(20) * xx[0] * xx[1] * xx[1] * xx[1]; };
std::function< real_t( const hyteg::Point3D& ) > exactV = []( const hyteg::Point3D& xx ) { return real_c(5) * xx[0] * xx[0] * xx[0] * xx[0] - real_c(5) * xx[1] * xx[1] * xx[1] * xx[1]; };
std::function< real_t( const hyteg::Point3D& ) > exactP = []( const hyteg::Point3D& xx ) { return real_c(60) * std::pow( xx[0], 2.0 ) * xx[1] - real_c(20) * std::pow( xx[1], 3.0 ); };
std::function< real_t( const hyteg::Point3D& ) > zero = []( const hyteg::Point3D& ) { return real_c(0); };
x.uvw.u.interpolate( exactU, level, hyteg::DirichletBoundary );
x.uvw.v.interpolate( exactV, level, hyteg::DirichletBoundary );
x_exact.uvw.u.interpolate( exactU, level );
x_exact.uvw.v.interpolate( exactV, level );
x_exact.p.interpolate( exactP, level );
// VTKOutput vtkOutput("../../output", "P2P1Stokes2DSchurCGConvergence", storage);
// vtkOutput.add( x.u );
// vtkOutput.add( x.v );
// vtkOutput.add( x.p );
// vtkOutput.add( x_exact.u );
// vtkOutput.add( x_exact.v );
// vtkOutput.add( x_exact.p );
// vtkOutput.add( err.u );
// vtkOutput.add( err.v );
// vtkOutput.add( err.p );
// vtkOutput.add( b.u );
// vtkOutput.add( b.v );
// vtkOutput.add( b.p );
// vtkOutput.write( level, 0 );
bool test2D = !storage->hasGlobalCells();
if ( test2D )
{
std::function< real_t( const hyteg::Point3D& ) > exactU = []( const hyteg::Point3D& xx ) {
return real_c( 20 ) * xx[0] * xx[1] * xx[1] * xx[1];
};
std::function< real_t( const hyteg::Point3D& ) > exactV = []( const hyteg::Point3D& xx ) {
return real_c( 5 ) * xx[0] * xx[0] * xx[0] * xx[0] - real_c( 5 ) * xx[1] * xx[1] * xx[1] * xx[1];
};
std::function< real_t( const hyteg::Point3D& ) > exactP = []( const hyteg::Point3D& xx ) {
return real_c( 60 ) * std::pow( xx[0], 2.0 ) * xx[1] - real_c( 20 ) * std::pow( xx[1], 3.0 );
};
x.uvw.interpolate( {exactU, exactV}, level, hyteg::DirichletBoundary );
x_exact.uvw.interpolate( {exactU, exactV}, level );
x_exact.p.interpolate( exactP, level );
}
else
{
std::function< real_t( const hyteg::Point3D& ) > exactU = []( const hyteg::Point3D& xx ) {
return -real_c( 4 ) * cos( real_c( 4 ) * xx[2] );
};
std::function< real_t( const hyteg::Point3D& ) > exactV = []( const hyteg::Point3D& xx ) {
return real_c( 8 ) * cos( real_c( 8 ) * xx[0] );
};
std::function< real_t( const hyteg::Point3D& ) > exactW = []( const hyteg::Point3D& xx ) {
return -real_c( 2 ) * cos( real_c( 2 ) * xx[1] );
};
x.uvw.interpolate( {exactU, exactV, exactW}, level, hyteg::DirichletBoundary );
x_exact.uvw.interpolate( {exactU, exactV, exactW}, level );
std::function< real_t( const hyteg::Point3D& ) > exactP = []( const hyteg::Point3D& xx ) {
return sin( real_c( 4 ) * xx[0] ) * sin( real_c( 8 ) * xx[1] ) * sin( real_c( 2 ) * xx[2] );
};
x_exact.p.interpolate( exactP, level );
std::function< real_t( const hyteg::Point3D& ) > rhs0 = []( const hyteg::Point3D& xx ) {
return real_c( 4 ) * cos( real_c( 4 ) * xx[0] ) * sin( real_c( 8 ) * xx[1] ) * sin( real_c( 2 ) * xx[2] ) - real_c(64) * cos( real_c(4) * xx[2] );
};
std::function< real_t( const hyteg::Point3D& ) > rhs1 = []( const hyteg::Point3D& xx ) {
return real_c( 8 ) * sin( real_c( 4 ) * xx[0] ) * cos( real_c( 8 ) * xx[1] ) * sin( real_c( 2 ) * xx[2] ) + real_c(512) * cos( real_c(8) * xx[0] );
};
std::function< real_t( const hyteg::Point3D& ) > rhs2 = []( const hyteg::Point3D& xx ) {
return real_c( 2 ) * sin( real_c( 4 ) * xx[0] ) * sin( real_c( 8 ) * xx[1] ) * cos( real_c( 2 ) * xx[2] ) - real_c(8) * cos( real_c(2) * xx[1] );
};
hyteg::P2VectorFunction< real_t > rhs( "rhs", storage, minLevel, level );
rhs.interpolate( {rhs0, rhs1, rhs2}, level );
P2ConstantMassOperator mass( storage, level, level );
for( uint_t k = 0; k < 3; k++ ) {
mass.apply( rhs[k], b.uvw[k], level, All );
}
}
uint_t localDoFs1 = hyteg::numberOfLocalDoFs< P2P1TaylorHoodFunctionTag >( *storage, level );
uint_t globalDoFs1 = hyteg::numberOfGlobalDoFs< P2P1TaylorHoodFunctionTag >( *storage, level );
......@@ -105,7 +134,8 @@ void P2P1SchurCGConvergenceTest( const uint_t & level, const MeshInfo & meshInfo
auto restriction = std::make_shared< P2toP2QuadraticRestriction >();
auto velocitySolver = std::make_shared< GeometricMultigridSolver< P2ConstantLaplaceOperator > >( storage, smoother, coarseGrid, restriction, prolongation, minLevel, level, 2, 2 );
auto loop = std::make_shared< SolverLoop< P2ConstantLaplaceOperator > >( velocitySolver, 10 );
StokesPCGSolverOld< P2P1TaylorHoodStokesOperator > solver( storage, loop, minLevel, level, 1e-10, 100, Inner | NeumannBoundary );
uint_t maxIter = test2D ? 100 : 2;
StokesPCGSolverOld< P2P1TaylorHoodStokesOperator > solver( storage, loop, minLevel, level, 1e-10, maxIter, Inner | NeumannBoundary );
walberla::WcTimer timer;
solver.solve( A, x, b, level );
......@@ -119,35 +149,57 @@ void P2P1SchurCGConvergenceTest( const uint_t & level, const MeshInfo & meshInfo
err.assign( {1.0, -1.0}, {x, x_exact}, level );
bool vtkOut = false;
if ( vtkOut )
{
VTKOutput vtkOutput( "../../output", "P2P1StokesPCGSolverConvergence", storage );
vtkOutput.add( x.uvw );
vtkOutput.add( x.p );
vtkOutput.add( x_exact.uvw );
vtkOutput.add( x_exact.p );
vtkOutput.add( err.uvw );
vtkOutput.add( err.p );
vtkOutput.add( b.uvw );
vtkOutput.add( b.p );
vtkOutput.write( level, 0 );
}
real_t discr_l2_err_u = std::sqrt( err.uvw.u.dotGlobal( err.uvw.u, level ) / (real_t) globalDoFsPerVelocityComponent );
real_t discr_l2_err_v = std::sqrt( err.uvw.v.dotGlobal( err.uvw.v, level ) / (real_t) globalDoFsPerVelocityComponent );
real_t discr_l2_err_p = std::sqrt( err.p.dotGlobal( err.p, level ) / (real_t) globalDoFsPressure );
real_t residuum_l2_1 = std::sqrt( residuum.dotGlobal( residuum, level, Inner ) / (real_t) globalDoFs1 );
WALBERLA_LOG_INFO_ON_ROOT( "discrete L2 error u = " << discr_l2_err_u );
WALBERLA_LOG_INFO_ON_ROOT( "discrete L2 error v = " << discr_l2_err_v );
if ( !test2D )
{
real_t discr_l2_err_w = std::sqrt( err.uvw.w.dotGlobal( err.uvw.w, level ) / (real_t) globalDoFsPerVelocityComponent );
WALBERLA_LOG_INFO_ON_ROOT( "discrete L2 error w = " << discr_l2_err_w );
}
WALBERLA_LOG_INFO_ON_ROOT( "discrete L2 error p = " << discr_l2_err_p );
WALBERLA_LOG_INFO_ON_ROOT( "residuum 1 = " << residuum_l2_1 );
// vtkOutput.write( level, 3 );
WALBERLA_CHECK_LESS( discr_l2_err_u, 3.5e-03 );
WALBERLA_CHECK_LESS( discr_l2_err_v, 2.4e-03 );
WALBERLA_CHECK_LESS( discr_l2_err_p, 2.9e-01 );
WALBERLA_CHECK_LESS( residuum_l2_1, 2.5e-09 );
// vtkOutput.write( level, 3 );
if ( test2D )
{
WALBERLA_CHECK_LESS( discr_l2_err_u, 3.5e-03 );
WALBERLA_CHECK_LESS( discr_l2_err_v, 2.4e-03 );
WALBERLA_CHECK_LESS( discr_l2_err_p, 2.9e-01 );
WALBERLA_CHECK_LESS( residuum_l2_1, 2.5e-09 );
}
}
}
using namespace hyteg;
int main( int argc, char* argv[] )
{
walberla::Environment walberlaEnv( argc, argv );
walberla::MPIManager::instance()->useWorldComm();
walberla::Environment walberlaEnv( argc, argv );
walberla::MPIManager::instance()->useWorldComm();
P2P1SchurCGConvergenceTest( 3, hyteg::MeshInfo::fromGmshFile( "../../data/meshes/quad_center_at_origin_4el.msh" ) );
P2P1PCGSolverConvergenceTest( 3, hyteg::MeshInfo::fromGmshFile( "../../data/meshes/quad_center_at_origin_4el.msh" ) );
P2P1PCGSolverConvergenceTest( 2, hyteg::MeshInfo::meshCuboid( Point3D( {0.0, 0.0, 0.0} ), Point3D( {1.0, 1.0, 1.0} ), 1, 1, 1 ) );
return EXIT_SUCCESS;
return EXIT_SUCCESS;
}
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