Commit 55e3ea1b authored by Nils Kohl's avatar Nils Kohl 🌝
Browse files

Fixes after rebasing pressure-cg branch (60d99b0e) onto master (6da4570a)

parent 15a5f5b2
Pipeline #27400 failed with stages
in 15 minutes and 26 seconds
......@@ -256,7 +256,7 @@ std::shared_ptr< SetupPrimitiveStorage >
};
setupStorage->setMeshBoundaryFlagsOnBoundary( 0, 0, true );
setupStorage->setMeshBoundaryFlagsByVertexLocation( 1, freeslip );
setupStorage->setMeshBoundaryFlagsByVertexLocation( 3, freeslip );
setupStorage->setMeshBoundaryFlagsByVertexLocation( 1, inflow );
AnnulusMap::setMap( *setupStorage );
......@@ -315,9 +315,9 @@ class SolverAsOperator : public Operator< typename OperatorType::dstType, typena
copySrc.assign( { 1 }, { src }, level, All );
copySrc.p.interpolate( 0, level, All );
solver_.solve( op_, copyDst, copySrc, level );
dst.u.assign( { 1 }, { copyDst.u }, level, All );
dst.v.assign( { 1 }, { copyDst.v }, level, All );
dst.w.assign( { 1 }, { copyDst.w }, level, All );
dst.uvw.u.assign( { 1 }, { copyDst.uvw.u }, level, All );
dst.uvw.v.assign( { 1 }, { copyDst.uvw.v }, level, All );
dst.uvw.w.assign( { 1 }, { copyDst.uvw.w }, level, All );
}
private:
......@@ -357,9 +357,11 @@ class StokesVelocityPressureBlockPreconditioner : public Solver< OperatorType >
};
template < typename StokesFunctionType,
typename StokesFunctionTag,
typename StokesOperatorType,
typename ProjectNormalOperatorType >
typename StokesFunctionTag,
typename StokesOperatorType,
typename ProjectNormalOperatorType,
typename RestrictionType,
typename ProlongationType >
void run( std::shared_ptr< walberla::config::Config > cfg )
{
const walberla::Config::BlockHandle mainConf = cfg->getBlock( "Parameters" );
......@@ -442,7 +444,6 @@ void run( std::shared_ptr< walberla::config::Config > cfg )
vtkOutput.add( u );
}
// using StokesOperator = hyteg::StrongFreeSlipWrapper< hyteg::P1StokesOperator, hyteg::P1ProjectNormalOperator >;
using StokesOperatorFS = hyteg::StrongFreeSlipWrapper< StokesOperatorType, ProjectNormalOperatorType >;
auto stokes = std::make_shared< StokesOperatorType >( storage, minLevel, maxLevel );
auto normalsRect = []( auto, Point3D& n ) { n = Point3D( {0, -1} ); };
......@@ -484,9 +485,8 @@ void run( std::shared_ptr< walberla::config::Config > cfg )
s->setPrintInfo( true );
using StokesFunction = P1StokesFunction< real_t >;
using StokesBlockOperator =
hyteg::block_structure::StokesFreeSlipBlockOperator< P1BlendingStokesOperator, hyteg::P1ProjectNormalOperator >;
hyteg::block_structure::StokesFreeSlipBlockOperator< StokesOperatorType, ProjectNormalOperatorType >;
auto stokesOp = std::make_shared< StokesBlockOperator >( stokes, projection, false );
// using StokesBlockOperator = hyteg::block_structure::StokesBlockOperator< P1BlendingStokesOperator >;
// auto stokesOp = std::make_shared< StokesBlockOperator >( stokes );
......@@ -502,7 +502,7 @@ void run( std::shared_ptr< walberla::config::Config > cfg )
}
{
P1StokesFunction< real_t > uu( "uu", storage, minLevel, maxLevel );
StokesFunctionType uu( "uu", storage, minLevel, maxLevel );
if ( useRectangleScenario )
applyDirichletBCRectangle( channelLength, channelHeight, maxLevel, uu );
......@@ -511,27 +511,23 @@ void run( std::shared_ptr< walberla::config::Config > cfg )
vtkOutput.add( uu );
using AType = StokesBlockOperator::AType;
using ProjectionType = P1ProjectNormalOperator;
using RestrictionType = P1P1StokesToP1P1StokesRestriction;
using ProlongationType = P1P1StokesToP1P1StokesProlongation;
using AType = typename StokesBlockOperator::AType;
// restriction:
// using FreeSlipRestrictionType = StrongFreeSlipRestrictionWrapper< StokesFunction, RestrictionType, ProjectionType >;
// auto restriction = std::make_shared< RestrictionType >();
// auto restrictionFS = std::make_shared< FreeSlipRestrictionType >( restriction, projection, FreeslipBoundary );
using FreeSlipRestrictionType = StrongFreeSlipRestrictionWrapper< StokesFunctionType, RestrictionType, ProjectNormalOperatorType >;
auto restriction = std::make_shared< RestrictionType >();
auto restrictionFS = std::make_shared< FreeSlipRestrictionType >( restriction, projection, FreeslipBoundary );
using FreeSlipRestrictionType = RestrictionType;
auto restrictionFS = std::make_shared< RestrictionType >();
// using FreeSlipRestrictionType = RestrictionType;
// auto restrictionFS = std::make_shared< RestrictionType >();
// prolongation:
// using FreeSlipProlongationType = StrongFreeSlipProlongationWrapper< StokesFunction, ProlongationType, ProjectionType >;
// auto prolongation = std::make_shared< ProlongationType >();
// auto prolongationFS = std::make_shared< FreeSlipProlongationType >( prolongation, projection, FreeslipBoundary );
using FreeSlipProlongationType = StrongFreeSlipProlongationWrapper< StokesFunctionType, ProlongationType, ProjectNormalOperatorType >;
auto prolongation = std::make_shared< ProlongationType >();
auto prolongationFS = std::make_shared< FreeSlipProlongationType >( prolongation, projection, FreeslipBoundary );
using FreeSlipProlongationType = ProlongationType;
auto prolongationFS = std::make_shared< ProlongationType >();
// using FreeSlipProlongationType = ProlongationType;
// auto prolongationFS = std::make_shared< ProlongationType >();
// coarse grid solver:
// using FreeSlipCoarseGridSolverType = StrongFreeSlipSolverWrapper< AType, CGSolver< AType >, ProjectionType >;
......@@ -568,7 +564,7 @@ void run( std::shared_ptr< walberla::config::Config > cfg )
auto mgSolverFS = std::make_shared< MGType > ( mgSolverFS_, projection );
*/
using MGType = StrongFreeSlipSolverWrapper< AType, Repeater< AType >, ProjectionType >;
using MGType = StrongFreeSlipSolverWrapper< AType, Repeater< AType >, ProjectNormalOperatorType >;
auto mgSolverInner = std::make_shared< GeometricMultigridSolver< AType > >(
storage, smootherFS, coarseGridSolverFS, restrictionFS, prolongationFS, minLevel, maxLevel, 2, 2 );
auto mgSolver = std::make_shared< Repeater< AType > >( stokesOp->getA(), mgSolverInner, 1e-11, 100 );
......@@ -639,15 +635,17 @@ int main( int argc, char* argv[] )
run< P1StokesFunction< real_t >, // function type
P1StokesFunction< real_t >::Tag, // tag
P1BlendingStokesOperator, // operator
P1ProjectNormalOperator // projection
>( cfg );
P1ProjectNormalOperator, // projection
P1P1StokesToP1P1StokesRestriction,
P1P1StokesToP1P1StokesProlongation >( cfg );
}
else
{
run< P2P1TaylorHoodFunction< real_t >, // function type
P2P1TaylorHoodFunction< real_t >::Tag, // tag
P2P1TaylorHoodStokesOperator, // operator
P2ProjectNormalOperator // projection
>( cfg );
P2ProjectNormalOperator, // projection
P2P1StokesToP2P1StokesRestriction,
P2P1StokesToP2P1StokesProlongation >( cfg );
}
}
......@@ -47,7 +47,7 @@ class P2P1TaylorHoodStokesOperator : public Operator< P2P1TaylorHoodFunction< re
, divT_y( storage, minLevel, maxLevel )
, divT_z( storage, minLevel, maxLevel )
, divT( storage, minLevel, maxLevel )
, pspg_( storage, minLevel, maxLevel )
, pspg( storage, minLevel, maxLevel )
, pspg_inv_diag_( storage, minLevel, maxLevel )
, hasGlobalCells_( storage->hasGlobalCells() )
{}
......
......@@ -207,21 +207,21 @@ class EpsilonBlockOperator : public Operator< typename OperatorType::srcType, ty
const DoFType flag,
const UpdateType updateType = Replace ) const
{
get_A_uu().apply( src.u, dst.u, level, flag, updateType );
get_A_uv().apply( src.v, dst.u, level, flag, Add );
get_A_uu().apply( src.uvw.u, dst.uvw.u, level, flag, updateType );
get_A_uv().apply( src.uvw.v, dst.uvw.u, level, flag, Add );
if ( hasGlobalCells_ )
get_A_uw().apply( src.w, dst.u, level, flag, Add );
get_A_uw().apply( src.uvw.w, dst.uvw.u, level, flag, Add );
get_A_vu().apply( src.u, dst.v, level, flag, updateType );
get_A_vv().apply( src.v, dst.v, level, flag, Add );
get_A_vu().apply( src.uvw.u, dst.uvw.v, level, flag, updateType );
get_A_vv().apply( src.uvw.v, dst.uvw.v, level, flag, Add );
if ( hasGlobalCells_ )
get_A_vw().apply( src.w, dst.v, level, flag, Add );
get_A_vw().apply( src.uvw.w, dst.uvw.v, level, flag, Add );
if ( hasGlobalCells_ )
{
get_A_wu().apply( src.u, dst.w, level, flag, updateType );
get_A_wv().apply( src.v, dst.w, level, flag, Add );
get_A_ww().apply( src.w, dst.w, level, flag, Add );
get_A_wu().apply( src.uvw.u, dst.uvw.w, level, flag, updateType );
get_A_wv().apply( src.uvw.v, dst.uvw.w, level, flag, Add );
get_A_ww().apply( src.uvw.w, dst.uvw.w, level, flag, Add );
}
}
......@@ -265,8 +265,8 @@ class DivTBlockOperator : public Operator< typename OperatorType::srcType, typen
void apply( const srcType& src, const dstType& dst, const uint_t level, const DoFType flag, const UpdateType updateType ) const
{
op_->divT_x.apply( src.p, dst.u, level, flag, updateType );
op_->divT_y.apply( src.p, dst.v, level, flag, updateType );
op_->divT_x.apply( src.p, dst.uvw.u, level, flag, updateType );
op_->divT_y.apply( src.p, dst.uvw.v, level, flag, updateType );
}
const std::shared_ptr< OperatorType >& op_;
......@@ -292,8 +292,8 @@ class DivBlockOperator : public Operator< typename OperatorType::srcType, typena
void apply( const srcType& src, const dstType& dst, const uint_t level, const DoFType flag, const UpdateType updateType ) const
{
op_->div_x.apply( src.u, dst.p, level, flag | DirichletBoundary, updateType );
op_->div_y.apply( src.v, dst.p, level, flag | DirichletBoundary, Add );
op_->div_x.apply( src.uvw.u, dst.p, level, flag | DirichletBoundary, updateType );
op_->div_y.apply( src.uvw.v, dst.p, level, flag | DirichletBoundary, Add );
}
const std::shared_ptr< OperatorType >& op_;
......
......@@ -53,24 +53,24 @@ class VectorialGaussSeidelSmoother : public Solver< OperatorType >
void solve( const OperatorType& A, const srcType& x, const dstType& b, const walberla::uint_t level ) override
{
A.get_A_uv().apply( x.v, tmp_.u, level, flag_, Replace );
A.get_A_uv().apply( x.uvw.v, tmp_.uvw.u, level, flag_, Replace );
if ( hasGlobalCells_ )
A.get_A_uw().apply( x.w, tmp_.u, level, flag_, Add );
tmp_.u.assign( { 1, -1 }, { b.u, tmp_.u }, level, flag_ );
A.get_A_uu().smooth_gs( x.u, tmp_.u, level, flag_ );
A.get_A_uw().apply( x.uvw.w, tmp_.uvw.u, level, flag_, Add );
tmp_.uvw.u.assign( { 1, -1 }, { b.uvw.u, tmp_.uvw.u }, level, flag_ );
A.get_A_uu().smooth_gs( x.uvw.u, tmp_.uvw.u, level, flag_ );
A.get_A_vu().apply( x.u, tmp_.v, level, flag_, Replace );
A.get_A_vu().apply( x.uvw.u, tmp_.uvw.v, level, flag_, Replace );
if ( hasGlobalCells_ )
A.get_A_vw().apply( x.w, tmp_.v, level, flag_, Add );
tmp_.v.assign( { 1, -1 }, { b.v, tmp_.v }, level, flag_ );
A.get_A_vv().smooth_gs( x.v, tmp_.v, level, flag_ );
A.get_A_vw().apply( x.uvw.w, tmp_.uvw.v, level, flag_, Add );
tmp_.uvw.v.assign( { 1, -1 }, { b.uvw.v, tmp_.uvw.v }, level, flag_ );
A.get_A_vv().smooth_gs( x.uvw.v, tmp_.uvw.v, level, flag_ );
if ( hasGlobalCells_ )
{
A.get_A_wu().apply( x.u, tmp_.w, level, flag_, Replace );
A.get_A_wv().apply( x.v, tmp_.w, level, flag_, Add );
tmp_.w.assign( { 1, -1 }, { b.w, tmp_.w }, level, flag_ );
A.get_A_ww().smooth_gs( x.w, tmp_.w, level, flag_ );
A.get_A_wu().apply( x.uvw.u, tmp_.uvw.w, level, flag_, Replace );
A.get_A_wv().apply( x.uvw.v, tmp_.uvw.w, level, flag_, Add );
tmp_.uvw.w.assign( { 1, -1 }, { b.uvw.w, tmp_.uvw.w }, level, flag_ );
A.get_A_ww().smooth_gs( x.uvw.w, tmp_.uvw.w, level, flag_ );
}
}
......@@ -107,32 +107,32 @@ class VectorialJacobiSmoother : public Solver< OperatorType >
void solve( const OperatorType& A, const srcType& x, const dstType& b, const walberla::uint_t level ) override
{
A.get_A_uv().apply( x.v, tmp_.u, level, flag_, Replace );
A.get_A_uv().apply( x.uvw.v, tmp_.uvw.u, level, flag_, Replace );
if ( hasGlobalCells_ )
A.get_A_uw().apply( x.w, tmp_.u, level, flag_, Add );
tmp_.u.assign( { 1, -1 }, { b.u, tmp_.u }, level, flag_ );
A.get_A_uw().apply( x.uvw.w, tmp_.uvw.u, level, flag_, Add );
tmp_.uvw.u.assign( { 1, -1 }, { b.u, tmp_.uvw.u }, level, flag_ );
A.get_A_vu().apply( x.u, tmp_.v, level, flag_, Replace );
A.get_A_vu().apply( x.uvw.u, tmp_.uvw.v, level, flag_, Replace );
if ( hasGlobalCells_ )
A.get_A_vw().apply( x.w, tmp_.v, level, flag_, Add );
tmp_.v.assign( { 1, -1 }, { b.v, tmp_.v }, level, flag_ );
A.get_A_vw().apply( x.uvw.w, tmp_.uvw.v, level, flag_, Add );
tmp_.uvw.v.assign( { 1, -1 }, { b.v, tmp_.uvw.v }, level, flag_ );
if ( hasGlobalCells_ )
{
A.get_A_wu().apply( x.u, tmp_.w, level, flag_, Replace );
A.get_A_wv().apply( x.v, tmp_.w, level, flag_, Add );
tmp_.w.assign( { 1, -1 }, { b.w, tmp_.w }, level, flag_ );
A.get_A_wu().apply( x.uvw.u, tmp_.uvw.w, level, flag_, Replace );
A.get_A_wv().apply( x.uvw.v, tmp_.uvw.w, level, flag_, Add );
tmp_.uvw.w.assign( { 1, -1 }, { b.w, tmp_.uvw.w }, level, flag_ );
}
tmp2_.u.assign( { 1.0 }, { x.u }, level, All );
tmp2_.v.assign( { 1.0 }, { x.v }, level, All );
tmp2_.uvw.u.assign( { 1.0 }, { x.uvw.u }, level, All );
tmp2_.uvw.v.assign( { 1.0 }, { x.uvw.v }, level, All );
if ( hasGlobalCells_ )
tmp2_.w.assign( { 1.0 }, { x.w }, level, All );
tmp2_.uvw.w.assign( { 1.0 }, { x.uvw.w }, level, All );
A.get_A_uu().smooth_jac( x.u, tmp_.u, tmp2_.u, relax_, level, flag_ );
A.get_A_vv().smooth_jac( x.v, tmp_.v, tmp2_.v, relax_, level, flag_ );
A.get_A_uu().smooth_jac( x.uvw.u, tmp_.uvw.u, tmp2_.uvw.u, relax_, level, flag_ );
A.get_A_vv().smooth_jac( x.uvw.v, tmp_.uvw.v, tmp2_.uvw.v, relax_, level, flag_ );
if ( hasGlobalCells_ )
A.get_A_ww().smooth_jac( x.w, tmp_.w, tmp2_.w, relax_, level, flag_ );
A.get_A_ww().smooth_jac( x.uvw.w, tmp_.uvw.w, tmp2_.uvw.w, relax_, level, flag_ );
}
private:
......
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