Skip to content
GitLab
Commit 79b12b5c authored by Andreas Wagner's avatar Andreas Wagner
Browse files

experimenting with the uzawa api on free-slip operators

parent 53b0eb97
Pipeline #25954 failed with stages
in 5 seconds
Hide whitespace changes
Inline Side-by-side
apps/stokesFreeSlip/StokesFreeSlip.cpp
View file @ 79b12b5c
......@@ -45,6 +45,54 @@ using walberla::real_t;
using walberla::uint_t;
using namespace hyteg;
using StokesOperator = hyteg::StrongFreeSlipOperatorWrapper< P1StokesOperator, P1ProjectNormalOperator >;
// using VelocityOperator = hyteg::StrongFreeSlipOperatorWrapper< P1ConstantLaplaceOperator, P1ProjectNormalOperator >;
using VelocityOperator = P1ConstantLaplaceOperator;
using StokesRestrictionOperator = hyteg::StrongFreeSlipRestrictionWrapper< P1StokesFunction< real_t >, P1P1StokesToP1P1StokesRestriction, P1ProjectNormalOperator >;
using StokesProlongationOperator = hyteg::StrongFreeSlipProlongationWrapper< P1StokesFunction< real_t >, P1P1StokesToP1P1StokesProlongation, P1ProjectNormalOperator >;
//using StokesGSSmoother = hyteg::StrongFreeSlipSolverWrapper< VelocityOperator, GaussSeidelSmoother< VelocityOperator > , P1ProjectNormalOperator >;
using StokesGSSmoother = GaussSeidelSmoother< VelocityOperator >;
std::shared_ptr< Solver< StokesOperator > >
stokesGMGUzawaSolverForP1P1( const std::shared_ptr< PrimitiveStorage >& storage,
const uint_t& minLevel,
const uint_t& maxLevel,
const uint_t& preSmoothingSteps,
const uint_t& postSmoothingSteps,
const real_t& uzawaSmootherOmega,
const std::shared_ptr< P1ProjectNormalOperator > projectNormalOperator )
{
auto pressurePreconditioner =
std::make_shared< StokesPressureBlockPreconditioner< StokesOperator, P1LumpedInvMassOperator > >(
storage, minLevel, minLevel );
auto pressurePreconditionedMinResSolver =
std::make_shared< MinResSolver< StokesOperator > >( storage, minLevel, minLevel, 1000, 1e-12, pressurePreconditioner );
auto bareRestriction = std::make_shared< P1P1StokesToP1P1StokesRestriction > ();
auto stokesRestriction = std::make_shared< StokesRestrictionOperator >( bareRestriction, projectNormalOperator, FreeslipBoundary );
auto bareProlongation = std::make_shared< P1P1StokesToP1P1StokesProlongation > ();
auto stokesProlongation = std::make_shared< StokesProlongationOperator >( bareProlongation, projectNormalOperator, FreeslipBoundary );
auto gaussSeidel = std::make_shared< GaussSeidelSmoother< VelocityOperator > >();
//auto bareGaussSeidel = std::make_shared< GaussSeidelSmoother< VelocityOperator > >();
//auto gaussSeidel = std::make_shared< StokesGSSmoother >( bareGaussSeidel, projectNormalOperator, FreeslipBoundary );
auto uzawaVelocityPreconditioner =
std::make_shared< StokesVelocityBlockBlockDiagonalPreconditioner< StokesOperator > >( storage, gaussSeidel );
auto uzawaSmoother = std::make_shared< UzawaSmoother< StokesOperator > >(
storage, uzawaVelocityPreconditioner, minLevel, maxLevel, uzawaSmootherOmega );
auto gmgSolver = std::make_shared< GeometricMultigridSolver< StokesOperator > >( storage,
uzawaSmoother,
pressurePreconditionedMinResSolver,
stokesRestriction,
stokesProlongation,
minLevel,
maxLevel,
preSmoothingSteps,
postSmoothingSteps,
2 );
return gmgSolver;
}
std::shared_ptr< SetupPrimitiveStorage >
setupStorageRectangle( const double channelLength, const double channelHeight, const uint_t ny )
{
......@@ -227,7 +275,6 @@ int main( int argc, char* argv[] )
vtkOutput.add( u );
}
using StokesOperator = hyteg::StrongFreeSlipWrapper< hyteg::P1StokesOperator, hyteg::P1ProjectNormalOperator >;
auto stokes = std::make_shared< hyteg::P1StokesOperator >( storage, minLevel, maxLevel );
auto normalsRect = []( auto, Point3D& n ) { n = Point3D( { 0, -1 } ); };
auto normalsAnn = [=]( Point3D p, Point3D& n ) { n = Point3D( { p[0] / rmax, p[1] / rmax } ); };
......@@ -244,8 +291,18 @@ int main( int argc, char* argv[] )
StokesOperator L( stokes, projection, FreeslipBoundary );
auto solver = hyteg::solvertemplates::stokesMinResSolver< StokesOperator >( storage, maxLevel, 1e-6, 1000 );
solver->solve( L, u, f, maxLevel );
bool useDefaultMinres = false;
if (useDefaultMinres)
{
auto solver = hyteg::solvertemplates::stokesMinResSolver< StokesOperator >( storage, maxLevel, 1e-12, 1000 );
solver->solve( L, u, f, maxLevel );
}
else
{
auto solver = stokesGMGUzawaSolverForP1P1( storage, minLevel, maxLevel, 2, 2, 1./3., projection );
solver->solve( L, u, f, maxLevel );
}
if ( mainConf.getParameter< bool >( "VTKOutput" ) )
{
......
apps/stokesFreeSlip/StokesFreeSlip.prm
View file @ 79b12b5c
Parameters
{
minLevel 2;
maxLevel 2;
minLevel 4;
maxLevel 4;
// which scenario should be loaded?
useRectangleScenario false;
......
src/hyteg/composites/P1StokesOperator.hpp
View file @ 79b12b5c
......@@ -26,69 +26,196 @@
namespace hyteg {
class P1StokesOperator : public Operator< P1StokesFunction< real_t >, P1StokesFunction< real_t > >
class P1VectorialLaplaceOperator : public Operator< P1StokesFunction< real_t >, P1StokesFunction< real_t > >
{
public:
typedef P1ConstantLaplaceOperator VelocityOperator_T;
typedef P1ConstantLaplaceOperator PressureOperator_T;
typedef P1StokesBlockPreconditioner BlockPreconditioner_T;
using ScalarAType = P1ConstantLaplaceOperator;
P1StokesOperator( const std::shared_ptr< PrimitiveStorage >& storage, size_t minLevel, size_t maxLevel )
: Operator( storage, minLevel, maxLevel )
, A( storage, minLevel, maxLevel )
, div_x( storage, minLevel, maxLevel )
, div_y( storage, minLevel, maxLevel )
, div_z( storage, minLevel, maxLevel )
, divT_x( storage, minLevel, maxLevel )
, divT_y( storage, minLevel, maxLevel )
, divT_z( storage, minLevel, maxLevel )
, pspg( storage, minLevel, maxLevel )
, pspg_inv_diag_( storage, minLevel, maxLevel )
, hasGlobalCells_( storage->hasGlobalCells() )
P1VectorialLaplaceOperator( const std::shared_ptr< PrimitiveStorage >& storage, size_t minLevel, size_t maxLevel )
: Operator( storage, minLevel, maxLevel )
, A( storage, minLevel, maxLevel )
, hasGlobalCells_( storage->hasGlobalCells() )
{}
void apply( const P1StokesFunction< real_t >& src,
const P1StokesFunction< real_t >& dst,
const size_t level,
DoFType flag ) const
DoFType flag,
UpdateType updateType = Replace) const
{
WALBERLA_ASSERT_NOT_IDENTICAL( std::addressof( src ), std::addressof( dst ) );
A.apply( src.u, dst.u, level, flag, Replace );
divT_x.apply( src.p, dst.u, level, flag, Add );
A.apply( src.v, dst.v, level, flag, Replace );
divT_y.apply( src.p, dst.v, level, flag, Add );
A.apply( src.u, dst.u, level, flag, updateType );
A.apply( src.v, dst.v, level, flag, updateType );
if ( hasGlobalCells_ )
if (hasGlobalCells_)
{
A.apply( src.w, dst.w, level, flag, Replace );
divT_z.apply( src.p, dst.w, level, flag, Add );
A.apply( src.w, dst.w, level, flag, updateType );
}
}
const ScalarAType & getScalarA() const { return A; }
div_x.apply( src.u, dst.p, level, flag, Replace );
private:
ScalarAType A;
bool hasGlobalCells_;
};
class P1DivOperator : public Operator< P1StokesFunction< real_t >, P1StokesFunction< real_t > >
{
public:
P1DivOperator( const std::shared_ptr< PrimitiveStorage >& storage, size_t minLevel, size_t maxLevel )
: Operator( storage, minLevel, maxLevel )
, div_x( storage, minLevel, maxLevel )
, div_y( storage, minLevel, maxLevel )
, div_z( storage, minLevel, maxLevel )
, hasGlobalCells_( storage->hasGlobalCells() )
{}
void apply( const P1StokesFunction< real_t >& src,
const P1StokesFunction< real_t >& dst,
const size_t level,
DoFType flag,
UpdateType updateType = Replace) const
{
WALBERLA_ASSERT_NOT_IDENTICAL( std::addressof( src ), std::addressof( dst ) );
div_x.apply( src.u, dst.p, level, flag, updateType );
div_y.apply( src.v, dst.p, level, flag, Add );
if ( hasGlobalCells_ )
{
div_z.apply( src.w, dst.p, level, flag, Add );
}
pspg.apply( src.p, dst.p, level, flag, Add );
}
P1ConstantLaplaceOperator A;
private:
P1DivxOperator div_x;
P1DivyOperator div_y;
P1DivzOperator div_z;
bool hasGlobalCells_;
};
class P1DivTOperator : public Operator< P1StokesFunction< real_t >, P1StokesFunction< real_t > >
{
public:
P1DivTOperator( const std::shared_ptr< PrimitiveStorage >& storage, size_t minLevel, size_t maxLevel )
: Operator( storage, minLevel, maxLevel )
, divT_x( storage, minLevel, maxLevel )
, divT_y( storage, minLevel, maxLevel )
, divT_z( storage, minLevel, maxLevel )
, hasGlobalCells_( storage->hasGlobalCells() )
{}
void apply( const P1StokesFunction< real_t >& src,
const P1StokesFunction< real_t >& dst,
const size_t level,
DoFType flag,
UpdateType updateType = Replace) const
{
WALBERLA_ASSERT_NOT_IDENTICAL( std::addressof( src ), std::addressof( dst ) );
divT_x.apply( src.p, dst.u, level, flag, updateType );
divT_y.apply( src.p, dst.v, level, flag, updateType );
if ( hasGlobalCells_ )
{
divT_z.apply( src.p, dst.w, level, flag, updateType );
}
}
private:
P1DivTxOperator divT_x;
P1DivTyOperator divT_y;
P1DivTzOperator divT_z;
bool hasGlobalCells_;
};
class P1StabOperator: public Operator< P1StokesFunction< real_t >, P1StokesFunction< real_t > >
{
public:
P1StabOperator( const std::shared_ptr< PrimitiveStorage >& storage, size_t minLevel, size_t maxLevel )
: Operator( storage, minLevel, maxLevel )
, pspg( storage, minLevel, maxLevel )
, pspg_inv_diag_( storage, minLevel, maxLevel )
, hasGlobalCells_( storage->hasGlobalCells() )
{}
void apply( const P1StokesFunction< real_t >& src,
const P1StokesFunction< real_t >& dst,
const size_t level,
DoFType flag,
UpdateType updateType = Replace) const
{
WALBERLA_ASSERT_NOT_IDENTICAL( std::addressof( src ), std::addressof( dst ) );
pspg.apply( src.p, dst.p, level, flag, updateType );
}
void applyInvDiag( const P1StokesFunction< real_t >& src,
const P1StokesFunction< real_t >& dst,
const size_t level,
DoFType flag,
UpdateType updateType = Replace) const
{
WALBERLA_ASSERT_NOT_IDENTICAL( std::addressof( src ), std::addressof( dst ) );
pspg_inv_diag_.apply( src.p, dst.p, level, flag, updateType );
}
private:
P1PSPGOperator pspg;
P1PSPGInvDiagOperator pspg_inv_diag_;
bool hasGlobalCells_;
};
class P1StokesOperator : public Operator< P1StokesFunction< real_t >, P1StokesFunction< real_t > >
{
public:
typedef P1ConstantLaplaceOperator VelocityOperator_T;
typedef P1ConstantLaplaceOperator PressureOperator_T;
typedef P1StokesBlockPreconditioner BlockPreconditioner_T;
using AType = P1VectorialLaplaceOperator;
using BType = P1DivOperator;
using BTType = P1DivTOperator;
using CType = P1StabOperator;
P1StokesOperator( const std::shared_ptr< PrimitiveStorage >& storage, size_t minLevel, size_t maxLevel )
: Operator( storage, minLevel, maxLevel )
, A( storage, minLevel, maxLevel )
, B( storage, minLevel, maxLevel )
, BT( storage, minLevel, maxLevel )
, C( storage, minLevel, maxLevel )
{}
void apply( const P1StokesFunction< real_t >& src,
const P1StokesFunction< real_t >& dst,
const size_t level,
DoFType flag ) const
{
WALBERLA_ASSERT_NOT_IDENTICAL( std::addressof( src ), std::addressof( dst ) );
A.apply( src, dst, level, flag, Replace );
BT.apply(src, dst, level, flag, Add );
B.apply(src, dst, level, flag, Replace );
C.apply( src, dst, level, flag, Add );
}
const AType & getA() const { return A; }
const BType & getB() const { return B; }
const BTType & getBT() const { return BT; }
const CType & getC() const { return C; }
private:
AType A;
BType B;
BTType BT;
CType C;
};
template <>
struct has_pspg_block< P1StokesOperator >
{
......
src/hyteg/composites/P2P1TaylorHoodStokesOperator.hpp
View file @ 79b12b5c
......@@ -73,6 +73,8 @@ class P2P1TaylorHoodStokesOperator : public Operator< P2P1TaylorHoodFunction< re
}
}
const VelocityOperator_T & getA() const { return A; };
P2ConstantLaplaceOperator A;
P2ToP1ConstantDivxOperator div_x;
P2ToP1ConstantDivyOperator div_y;
......
src/hyteg/composites/StrongFreeSlipWrapper.hpp
View file @ 79b12b5c
......@@ -19,9 +19,14 @@
*/
#pragma once
// TODO: Move this into its own file
#include <type_traits>
#include "hyteg/Operator.hpp"
#include "hyteg/gridtransferoperators/ProlongationOperator.hpp"
#include "hyteg/gridtransferoperators/RestrictionOperator.hpp"
#include "hyteg/solvers/Solver.hpp"
namespace hyteg {
......@@ -43,10 +48,22 @@ namespace hyteg {
/// \tparam ProjOpType Type of the normal projection operator, which enforces free slip
/// \tparam PreProjection Should the projection operator also already be applied to the src vector? This is not necessary if it is already consistent!
template < typename OpType, typename ProjOpType, bool PreProjection = false >
class StrongFreeSlipWrapper : public Operator< typename OpType::srcType, typename OpType::dstType >
class StrongFreeSlipOperatorWrapper : public Operator< typename OpType::srcType, typename OpType::dstType >
{
public:
StrongFreeSlipWrapper( std::shared_ptr< OpType > op, std::shared_ptr< ProjOpType > projOp, DoFType projFlag )
using srcType = typename OpType::srcType;
using dstType = typename OpType::dstType;
typedef typename OpType::VelocityOperator_T VelocityOperator_T;
typedef typename OpType::PressureOperator_T PressureOperator_T;
typedef typename OpType::BlockPreconditioner_T BlockPreconditioner_T;
using AType = typename OpType::AType;
using BType = typename OpType::BType;
using BTType = typename OpType::BTType;
using CType = typename OpType::CType;
StrongFreeSlipOperatorWrapper( std::shared_ptr< OpType > op, std::shared_ptr< ProjOpType > projOp, DoFType projFlag )
: Operator< typename OpType::srcType, typename OpType::dstType >( op->getStorage(), op->getMinLevel(), op->getMaxLevel() )
, op_( op )
, projOp_( projOp )
......@@ -78,6 +95,17 @@ class StrongFreeSlipWrapper : public Operator< typename OpType::srcType, typenam
projOp_->apply( dst, level, projFlag_ );
}
void smooth_gs( const srcType & dst, const dstType & rhs, size_t level, DoFType flag ) const
{
op_->smooth_gs(dst, rhs, level, flag);
projOp_->apply( dst, level, projFlag_ );
}
const AType & getA() const { return op_->getA(); };
const BType & getB() const { return op_->getB(); };
const BTType & getBT() const { return op_->getBT(); };
const CType & getC() const { return op_->getC(); };
private:
std::shared_ptr< OpType > op_;
std::shared_ptr< ProjOpType > projOp_;
......@@ -87,4 +115,155 @@ class StrongFreeSlipWrapper : public Operator< typename OpType::srcType, typenam
std::shared_ptr< typename OpType::dstType > tmp_;
};
template < typename OpType, typename ProjOpType, bool PreProjection = false >
class StrongFreeSlipCompositeOperatorWrapper : public Operator< typename OpType::srcType, typename OpType::dstType >
{
public:
using srcType = typename OpType::srcType;
using dstType = typename OpType::dstType;
using AType = typename OpType::AType;
using BType = typename OpType::BType;
using BTType = typename OpType::BTType;
using CType = typename OpType::CType;
StrongFreeSlipCompositeOperatorWrapper( std::shared_ptr< OpType > op, std::shared_ptr< ProjOpType > projOp, DoFType projFlag )
: Operator< typename OpType::srcType, typename OpType::dstType >( op->getStorage(), op->getMinLevel(), op->getMaxLevel() )
, op_( op )
, projOp_( projOp )
, projFlag_( projFlag )
, tmp_( PreProjection ?
std::make_shared< typename OpType::srcType >( "tmp", op->getStorage(), op->getMinLevel(), op->getMaxLevel() ) :
nullptr )
{}
void apply( const typename OpType::srcType& src,
const typename OpType::dstType& dst,
size_t level,
DoFType flag,
UpdateType updateType = Replace ) const
{
WALBERLA_CHECK( updateType == Replace, "Operator concatenation only supported for updateType Replace" );
if ( PreProjection )
{
tmp_->assign( { 1 }, { src }, level, All );
projOp_->apply( *tmp_, level, projFlag_ );
op_->apply( *tmp_, dst, level, flag );
}
else
{
op_->apply( src, dst, level, flag );
}
projOp_->apply( dst, level, projFlag_ );
}
void smooth_gs( const srcType & dst, const dstType & rhs, size_t level, DoFType flag ) const
{
op_->smooth_gs(dst, rhs, level, flag);
projOp_->apply( dst, level, projFlag_ );
}
const AType & getA() const { return op_->getA(); };
const BType & getB() const { return op_->getB(); };
const BTType & getBT() const { return op_->getBT(); };
const CType & getC() const { return op_->getC(); };
private:
std::shared_ptr< OpType > op_;
std::shared_ptr< ProjOpType > projOp_;
DoFType projFlag_;
std::shared_ptr< typename OpType::dstType > tmp_;
};
// TODO: Avoid 1st argument by definition in projection operator!
template < typename FunctionType, typename ProlonType, typename ProjProlonType >
class StrongFreeSlipProlongationWrapper : public ProlongationOperator< FunctionType >
{
public:
StrongFreeSlipProlongationWrapper( std::shared_ptr< ProlonType > prolon,
std::shared_ptr< ProjProlonType > projOp,
DoFType projFlag )
: prolon_( prolon )
, projOp_( projOp )
, projFlag_( projFlag )
{}
void prolongate( const FunctionType& function, const uint_t& sourceLevel, const DoFType& flag ) const override
{
const uint_t destinationLevel = sourceLevel + 1;
prolon_->prolongate( function, sourceLevel, flag );
projOp_->apply( function, destinationLevel, projFlag_ );
}
private:
std::shared_ptr< ProlonType > prolon_;
std::shared_ptr< ProjProlonType > projOp_;
DoFType projFlag_;
std::shared_ptr< FunctionType > tmp_;
};
// TODO: Avoid 1st argument by definition in projection operator!
template < typename FunctionType, typename RestrType, typename ProjOpType >
class StrongFreeSlipRestrictionWrapper : public RestrictionOperator< FunctionType >
{
public:
StrongFreeSlipRestrictionWrapper( std::shared_ptr< RestrType > restr, std::shared_ptr< ProjOpType > projOp, DoFType projFlag )
: restr_( restr )
, projOp_( projOp )
, projFlag_( projFlag )
{}
void restrict( const FunctionType& function, const uint_t& sourceLevel, const DoFType& flag ) const override
{
const uint_t destinationLevel = sourceLevel - 1;
restr_->restrict( function, sourceLevel, flag );
projOp_->apply( function, destinationLevel, projFlag_ );
}
private:
std::shared_ptr< RestrType > restr_;
std::shared_ptr< ProjOpType > projOp_;
DoFType projFlag_;
std::shared_ptr< FunctionType > tmp_;
};
/// Enforces a compatible right-hand-side for an operator with a strong free-slip boundary condition.
template < typename OperatorType, typename SolverType, typename ProjSolverType >
class StrongFreeSlipSolverWrapper : public Solver< SolverType >
{
public:
using FunctionType = typename OperatorType::srcType;
StrongFreeSlipSolverWrapper( std::shared_ptr< SolverType > solver, std::shared_ptr< ProjSolverType > projOp, DoFType projFlag )
: solver_( solver )
, projOp_( projOp )
, projFlag_( projFlag )
{}
void solve( const OperatorType& op, const FunctionType& x, const FunctionType& b, uint_t level ) const override
{
// tmp_->assign( { 1 }, { b }, level, All );
// projOp_->apply( *tmp_, level, projFlag_ );
// solver_->solve( op, x, *tmp_, level );
}
private:
std::shared_ptr< SolverType > solver_;
std::shared_ptr< ProjSolverType > projOp_;
DoFType projFlag_;
std::shared_ptr< FunctionType > tmp_;
};
} // namespace hyteg
src/hyteg/solvers/UzawaSmoother.hpp
View file @ 79b12b5c
......@@ -133,15 +133,13 @@ class UzawaSmoother : public Solver< OperatorType >
std::false_type /* tensor */,
std::true_type /* PSPG */ ) const
{
A.divT_x.apply( x.p, r_.u, level, flag_, Replace );
r_.u.assign( {1.0, -1.0}, {b.u, r_.u}, level, flag_ );
A.getBT().apply(x, r_, level, flag_, Replace);
A.divT_y.apply( x.p, r_.v, level, flag_, Replace );
r_.u.assign( {1.0, -1.0}, {b.u, r_.u}, level, flag_ );
r_.v.assign( {1.0, -1.0}, {b.v, r_.v}, level, flag_ );
if ( hasGlobalCells_ )
{
A.divT_z.apply( x.p, r_.w, level, flag_, Replace );
r_.w.assign( {1.0, -1.0}, {b.w, r_.w}, level, flag_ );
}
......@@ -150,46 +148,40 @@ class UzawaSmoother : public Solver< OperatorType >
velocitySmoother_->solve( A, x, r_, level );
}
A.pspg.apply( x.p, r_.p, level, flag_, Replace );
A.div_x.apply( x.u, r_.p, level, flag_, Add );
A.div_y.apply( x.v, r_.p, level, flag_, Add );
if ( hasGlobalCells_ )
{
A.div_z.apply( x.w, r_.p, level, flag_, Add );
}
A.getC().apply( x, r_, level, flag_, Replace );
A.getB().apply( x, r_, level, flag_, Add );
r_.p.assign( {1.0, -1.0}, {b.p, r_.p}, level, flag_ );
r_.p.assign( {relaxParam_}, {r_.p}, level, flag_ );
A.pspg_inv_diag_.apply( r_.p, x.p, level, flag_, Add );
A.getC().applyInvDiag( r_, x, level, flag_, Add );