Add parent classes for smoothers
In our meeting today the question was posed, which methods we should include in the abstract Operator interface.
One idea was to include all the smoothing procedures to make it easier to switch between operators and smoothers at runtime. This is currently not possible and since all our operators have different smoothers, changing the operator often requires commenting in large amounts of code to pass compilation, even if the respective code paths are never taken.
This idea was rejected to avoid too much clutter in the base class in the long run and since from an idiomatic point of view a class should be able to implement all the methods of the parent class.
An alternative suggestion was to implement an interface for every smoother type. For the Jacobi, GS and symmetric GS this could look as follows:
template < typename Function >
class JacobiSmoothable
{
public:
virtual ~JacobiSmoothable() = default;
virtual void smooth_jac( const Function& dst,
const Function& rhs,
const Function& tmp,
const real_t& relax,
size_t level,
DoFType flag ) const = 0;
};
template < typename Function >
class GSSmoothable
{
public:
virtual ~GSSmoothable() = default;
virtual void smooth_gs( const Function& dst, const Function& rhs, size_t level, DoFType flag ) const = 0;
};
template < typename Function >
class GSBackwardsSmoothable
{
public:
virtual ~GSBackwardsSmoothable() = default;
virtual void smooth_gs_backwards( const Function& dst, const Function& rhs, size_t level, DoFType flag ) const = 0;
};
// ... all the other smoothers ..
Operators, which want to support some of the smoothers can then inherit from the respective interfaces. For a reduced version of the P1ConstantOperator which only support Jacobi and forward GS this could look like this:
template < class P1Form, bool Diagonal = false, bool Lumped = false, bool InvertDiagonal = false >
class P1ConstantOperator : public Operator< P1Function< real_t >, P1Function< real_t > >,
public JacobiSmoothable< P1Function< real_t > >,
public GSSmoothable< P1Function< real_t > >
{
public:
void smooth_jac( const P1Function< real_t >& dst,
const P1Function< real_t >& rhs,
const P1Function< real_t >& tmp,
const real_t& relax,
size_t level,
DoFType flag ) const override
{
// implementation ...
}
void smooth_gs( const P1Function< real_t >& dst, const P1Function< real_t >& rhs, size_t level, DoFType flag ) const override
{
// implementation ...
}
};
Our solvers can then dynamically at runtime check if the operator supports the given smoother, and throw an exception if this is not the case:
template < class OperatorType >
class GaussSeidelSmoother : public Solver< OperatorType >
{
public:
void solve( const OperatorType& A,
const typename OperatorType::srcType& x,
const typename OperatorType::dstType& b,
const size_t level ) override
{
if ( const auto* A_gs = dynamic_cast< const GSSmoothable< typename OperatorType::srcType >* >( &A ) )
{
A_gs->smooth_gs( x, b, level, flag_ );
}
else
{
throw std::runtime_error( "The Gauss-Seidel Operator requires the GSSmoothable interface." );
}
}
// ...
};
template < class OperatorType >
class SymmetricGaussSeidelSmoother : public Solver< OperatorType >
{
public:
void solve( const OperatorType& A,
const typename OperatorType::srcType& x,
const typename OperatorType::dstType& b,
const size_t level ) override
{
if ( const auto* A_gs = dynamic_cast< const GSSmoothable< typename OperatorType::srcType >* >( &A ) )
{
A_gs->smooth_gs( x, b, level, flag_ );
}
else
{
throw std::runtime_error( "The symmetric Gauss-Seidel Operator requires the GSSmoothable interface." );
}
if( const auto* A_gs_bw = dynamic_cast< const GSBackwardsSmoothable< typename OperatorType::srcType >* >( &A ))
{
A_gs_bw->smooth_gs_backwards( x, b, level, flag_ );
}
else
{
throw std::runtime_error( "The symmetric Gauss-Seidel Operator requires the GSBackwardsSmoothable interface." );
}
}
// ...
};
Thus in the example above the GaussSeidelSmoother
would work for the P1ConstantLaplaceOperator
, while the SymmetricGaussSeidelSmoother
would throw an exception:
P1ConstantLaplaceOperator A( /* ... */ );
GaussSeidelSmoother< P1ConstantLaplaceOperator > GS ( /* ... */ );
SymmetricGaussSeidelSmoother< P1ConstantLaplaceOperator > SGS ( /* ... */ );
// works
GS.solve( A, x, b, level );
std::cout << "GS finished!" << std::endl;
// throws an exception at runtime
SGS.solve( A, x, b, level );
std::cout << "SGS finished!" << std::endl;