Skip to content
GitLab
Commit 5bfaa2db authored by Andreas Wagner's avatar Andreas Wagner
Browse files

changes the HDGP1 operator interface to match the P1 operator interface

parent 5ddbd00d
Pipeline #25637 failed with stages
in 12 minutes and 58 seconds
Hide whitespace changes
Inline Side-by-side
apps/hdgp1/HDGP1Adaptivity.cpp
View file @ 5bfaa2db
......@@ -133,7 +133,7 @@ void runGMG( const walberla::shared_ptr< walberla::config::Config >& cfg )
laplaceOperator.getMacroEdgeForm().setKappa( kappa );
laplaceOperator.getMacroEdgeForm().setSigma( sigma );
hdgp1::setupScenarios::setupEpsOperator( dGScheme, laplaceOperator );
laplaceOperator.computeDiagonalOperatorValues( true );
laplaceOperator.computeInverseDiagonalOperatorValues();
hdgp1::HDGP1Function< real_t > smootherTmp1( "smoother tmp 1", storage, minLevel, maxLevel, localLevelInfo );
smootherTmp1.interpolate( 0, minLevel, DoFType::DirichletBoundary );
......@@ -141,9 +141,16 @@ void runGMG( const walberla::shared_ptr< walberla::config::Config >& cfg )
minLevel,
DoFType::Inner | DoFType::NeumannBoundary );
hdgp1::HDGP1Function< real_t > smootherTmp2( "smoother tmp 2", storage, minLevel, maxLevel, localLevelInfo );
hdgp1::HDGP1Function< real_t > smootherTmp3( "smoother tmp 3", storage, minLevel, maxLevel, localLevelInfo );
std::string smootherName;
auto smoother = hdgp1::setupScenarios::setupSmoother(
parameters, commMode, laplaceOperator, smootherName, *laplaceOperator.diagonalValues_, smootherTmp1, smootherTmp2 );
auto smoother = hdgp1::setupScenarios::setupSmoother( parameters,
commMode,
laplaceOperator,
smootherName,
*laplaceOperator.getInverseDiagonalValues(),
smootherTmp1,
smootherTmp2,
smootherTmp3 );
WALBERLA_LOG_INFO_ON_ROOT( "created smoother " << smootherName );
const auto pc = hdgp1::setupScenarios::setupCoarseGridPreconditioner( dGScheme );
......
apps/hdgp1/HDGP1AsymptoticErrorRate.cpp
View file @ 5bfaa2db
......@@ -112,9 +112,10 @@ void runGMG( const walberla::shared_ptr< walberla::config::Config >& cfg )
smTmp1.interpolate( 0, minLevel, DirichletBoundary );
smTmp1.interpolate( []( const Point3D& p ) { return p[0] * p[0] + p[1] * p[1]; }, minLevel, Inner | NeumannBoundary );
P1Function< real_t > smTmp2( "smoother tmp 2", storage, minLevel, maxLevel );
P1Function< real_t > smTmp3( "smoother tmp 3", storage, minLevel, maxLevel );
const auto smParameters = parameters.getOneBlock( "P1" );
auto smoother = hyteg::hdgp1::setupScenarios::setupSmoother(
smParameters, commMode, laplaceOperator, smootherNameP1, *laplaceOperator.diagonalValues_, smTmp1, smTmp2 );
smParameters, commMode, laplaceOperator, smootherNameP1, *laplaceOperator.getInverseDiagonalValues(), smTmp1, smTmp2, smTmp3 );
auto coarseGridSolver = std::make_shared< CGSolver< P1ElementwiseLaplaceOperator > >(
storage, minLevel, minLevel, maxCoarseIter, coarseTolerance );
......@@ -199,16 +200,17 @@ void runGMG( const walberla::shared_ptr< walberla::config::Config >& cfg )
WALBERLA_LOG_INFO_ON_ROOT( "starting hdG with " << smoothingSteps << " smoothing steps and sigma = " << sigma );
laplaceOperator.getMacroEdgeForm().setSigma( sigma );
laplaceOperator.computeDiagonalOperatorValues( true );
laplaceOperator.computeInverseDiagonalOperatorValues( );
hdgp1::HDGP1Function< real_t > smTmp1( "smoother tmp 1", storage, minLevel, maxLevel, localLevelInfo );
smTmp1.interpolate( 0, minLevel, DoFType::Inner );
smTmp1.interpolate(
[]( const Point3D& p, const Point3D& ) { return p[0] * p[0] + p[1] * p[1]; }, minLevel, DoFType::Inner );
hdgp1::HDGP1Function< real_t > smTmp2( "smoother tmp 1", storage, minLevel, maxLevel, localLevelInfo );
hdgp1::HDGP1Function< real_t > smTmp2( "smoother tmp 2", storage, minLevel, maxLevel, localLevelInfo );
hdgp1::HDGP1Function< real_t > smTmp3( "smoother tmp 3", storage, minLevel, maxLevel, localLevelInfo );
const auto smParameters = parameters.getOneBlock( "hdGP1" );
auto smoother = hdgp1::setupScenarios::setupSmoother(
smParameters, commMode, laplaceOperator, smootherNameHDG, *laplaceOperator.diagonalValues_, smTmp1, smTmp2 );
smParameters, commMode, laplaceOperator, smootherNameHDG, *laplaceOperator.getInverseDiagonalValues(), smTmp1, smTmp2, smTmp3 );
auto mgSolver = GeometricMultigridSolver< hdgp1::HDGP1LaplaceOperator >( storage,
tmpGeomFct,
......
apps/hdgp1/HDGP1Eigenmodes.cpp
View file @ 5bfaa2db
......@@ -73,7 +73,7 @@ void run( const walberla::shared_ptr< walberla::config::Config >& cfg )
const auto randFunc = [&]( const Point3D&, const Point3D& ) { return walberla::math::realRandom< real_t >(); };
hdgp1::HDGP1LaplaceOperator laplace( storage, level, level, localLevelInfo );
laplace.computeDiagonalOperatorValues( true );
laplace.computeInverseDiagonalOperatorValues( );
laplace.useMaxDiameterAsPenaltyEdgeLength();
laplace.getMacroEdgeForm().setSigma( sigma );
......@@ -84,6 +84,8 @@ void run( const walberla::shared_ptr< walberla::config::Config >& cfg )
WALBERLA_LOG_INFO( "sigma = " << sigma );
hdgp1::HDGP1Function< real_t > tmp( "tmp", storage, level, level, localLevelInfo );
real_t maxEigenvalue = -1;
if ( useMass )
{
......@@ -97,7 +99,7 @@ void run( const walberla::shared_ptr< walberla::config::Config >& cfg )
else
{
maxEigenvalue =
estimateSpectralRadiusWithPowerIteration( laplace, function, numPowerIter, storage, level, flag, localLevelInfo, flag );
estimateSpectralRadiusWithPowerIteration( laplace, function, tmp, numPowerIter, storage, level, flag );
}
WALBERLA_LOG_INFO( "maxEigenvalue = " << maxEigenvalue );
......
apps/hdgp1/HDGP1Eigenvalues.cpp
View file @ 5bfaa2db
......@@ -73,7 +73,9 @@ void runGMG( const walberla::shared_ptr< walberla::config::Config >& cfg )
function.interpolate( bc3, level, DoFType::DirichletBoundary );
hdgp1::HDGP1LaplaceOperator laplace( storage, level, level, localLevelInfo );
laplace.computeDiagonalOperatorValues( true );
laplace.computeInverseDiagonalOperatorValues();
hdgp1::HDGP1Function< real_t > tmp( "tmp", storage, level, level, localLevelInfo );
std::vector< real_t > sigmaList;
std::vector< real_t > spectralRadiusRicList;
......@@ -112,7 +114,7 @@ void runGMG( const walberla::shared_ptr< walberla::config::Config >& cfg )
[]( const Point3D& p, const Point3D& ) { return p[0] * p[0] + std::sin( 2 * pi * p[0] ); }, level, flag );
const auto spectralRadiusRic = estimateSpectralRadiusWithPowerIteration(
laplace, function, numPowerIterRic, storage, level, flag, localLevelInfo, flag );
laplace, function, tmp, numPowerIterRic, storage, level, flag );
spectralRadiusRicList.push_back( spectralRadiusRic );
WALBERLA_LOG_INFO( "richardson: sigma = " << sigma << " r = " << spectralRadiusRic );
......@@ -122,10 +124,10 @@ void runGMG( const walberla::shared_ptr< walberla::config::Config >& cfg )
function.interpolate( 0, level, DoFType::DirichletBoundary );
hdgp1::utils::fillRandomly( function, level, flag );
laplace.computeDiagonalOperatorValues( true );
const auto wrapper = hdgp1::opalg::multiply( *laplace.diagonalValues_, laplace );
laplace.computeInverseDiagonalOperatorValues( );
const auto wrapper = hdgp1::opalg::multiply( *laplace.getInverseDiagonalValues(), laplace );
const auto spectralRadiusJac = estimateSpectralRadiusWithPowerIteration(
wrapper, function, numPowerIterJac, storage, level, flag, localLevelInfo, flag );
wrapper, function, tmp, numPowerIterJac, storage, level, flag );
spectralRadiusJacList.push_back( spectralRadiusJac );
WALBERLA_LOG_INFO( "jacobi: sigma = " << sigma << " r = " << spectralRadiusJac );
......
apps/hdgp1/HDGP1GMG.cpp
View file @ 5bfaa2db
......@@ -167,17 +167,18 @@ void runGMG( const walberla::shared_ptr< walberla::config::Config >& cfg )
laplaceOperator.getFaceForm().setKappa( kappa );
laplaceOperator.getMacroEdgeForm().setSigma( sigma );
hdgp1::setupScenarios::setupEpsOperator( dGScheme, laplaceOperator );
laplaceOperator.computeDiagonalOperatorValues( true );
laplaceOperator.computeInverseDiagonalOperatorValues();
hdgp1::HDGP1Function< real_t > smootherTmp1( "smoother tmp 1", storage, minLevel, maxLevel, localLevelInfo );
smootherTmp1.interpolate( 0, minLevel, DoFType::DirichletBoundary );
smootherTmp1.interpolate( []( const Point3D& p, const Point3D& ) { return p[0] * p[0] + p[1] * p[1]; },
hdgp1::HDGP1Function< real_t > smTmp1( "smoother tmp 1", storage, minLevel, maxLevel, localLevelInfo );
smTmp1.interpolate( 0, minLevel, DoFType::DirichletBoundary );
smTmp1.interpolate( []( const Point3D& p, const Point3D& ) { return p[0] * p[0] + p[1] * p[1]; },
minLevel,
DoFType::Inner | DoFType::NeumannBoundary );
hdgp1::HDGP1Function< real_t > smootherTmp2( "smoother tmp 2", storage, minLevel, maxLevel, localLevelInfo );
hdgp1::HDGP1Function< real_t > smTmp2( "smoother tmp 2", storage, minLevel, maxLevel, localLevelInfo );
hdgp1::HDGP1Function< real_t > smTmp3( "smoother tmp 3", storage, minLevel, maxLevel, localLevelInfo );
std::string smootherName;
auto smoother = hdgp1::setupScenarios::setupSmoother(
parameters, commMode, laplaceOperator, smootherName, *laplaceOperator.diagonalValues_, smootherTmp1, smootherTmp2 );
parameters, commMode, laplaceOperator, smootherName, *laplaceOperator.getInverseDiagonalValues(), smTmp1, smTmp2, smTmp3 );
WALBERLA_LOG_INFO_ON_ROOT( "created smoother " << smootherName );
const auto pc = hdgp1::setupScenarios::setupCoarseGridPreconditioner( dGScheme );
......
apps/hdgp1/HDGP1TwoGrids.cpp
View file @ 5bfaa2db
......@@ -106,17 +106,18 @@ void runGMG( const walberla::shared_ptr< walberla::config::Config >& cfg )
laplaceOperator.getMacroEdgeForm().setKappa( kappa );
laplaceOperator.getMacroEdgeForm().setSigma( sigma );
hdgp1::setupScenarios::setupEpsOperator( dGScheme, laplaceOperator );
laplaceOperator.computeDiagonalOperatorValues( true );
laplaceOperator.computeInverseDiagonalOperatorValues();
hdgp1::HDGP1Function< real_t > smootherTmp1( "smoother tmp 1", storage, minLevel, maxLevel, localLevelInfo );
smootherTmp1.interpolate( 0, minLevel, DoFType::DirichletBoundary );
smootherTmp1.interpolate( []( const Point3D& p, const Point3D& ) { return p[0] * p[0] + p[1] * p[1]; },
hdgp1::HDGP1Function< real_t > smTmp1( "smoother tmp 1", storage, minLevel, maxLevel, localLevelInfo );
smTmp1.interpolate( 0, minLevel, DoFType::DirichletBoundary );
smTmp1.interpolate( []( const Point3D& p, const Point3D& ) { return p[0] * p[0] + p[1] * p[1]; },
minLevel,
DoFType::Inner | DoFType::NeumannBoundary );
hdgp1::HDGP1Function< real_t > smootherTmp2( "smoother tmp 2", storage, minLevel, maxLevel, localLevelInfo );
hdgp1::HDGP1Function< real_t > smTmp2( "smoother tmp 2", storage, minLevel, maxLevel, localLevelInfo );
hdgp1::HDGP1Function< real_t > smTmp3( "smoother tmp 3", storage, minLevel, maxLevel, localLevelInfo );
std::string smootherName;
auto smoother = hdgp1::setupScenarios::setupSmoother(
parameters, commMode, laplaceOperator, smootherName, *laplaceOperator.diagonalValues_, smootherTmp1, smootherTmp2 );
parameters, commMode, laplaceOperator, smootherName, *laplaceOperator.getInverseDiagonalValues(), smTmp1, smTmp2, smTmp3 );
WALBERLA_LOG_INFO_ON_ROOT( "created smoother " << smootherName );
const auto pc = hdgp1::setupScenarios::setupCoarseGridPreconditioner( dGScheme );
......
apps/hdgp1/P1GMGCom.cpp
View file @ 5bfaa2db
......@@ -98,14 +98,15 @@ int main( int argc, char** argv )
massOperator.apply( tmp, rightHandSide, maxLevel, DoFType::All, UpdateType::Replace );
// create smoothers
P1Function< real_t > smootherTmp1( "smoother tmp 1", storage, minLevel, maxLevel );
smootherTmp1.interpolate( 0, minLevel, DoFType::DirichletBoundary );
smootherTmp1.interpolate(
P1Function< real_t > smTmp1( "smoother tmp 1", storage, minLevel, maxLevel );
smTmp1.interpolate( 0, minLevel, DoFType::DirichletBoundary );
smTmp1.interpolate(
[]( const Point3D& p ) { return p[0] * p[0] + p[1] * p[1]; }, minLevel, DoFType::Inner | DoFType::NeumannBoundary );
P1Function< real_t > smootherTmp2( "smoother tmp 2", storage, minLevel, maxLevel );
P1Function< real_t > smTmp2( "smoother tmp 2", storage, minLevel, maxLevel );
P1Function< real_t > smTmp3( "smoother tmp 3", storage, minLevel, maxLevel );
std::string smootherName;
auto smoother = hdgp1::setupScenarios::setupSmoother(
parameters, comm, laplaceOperator, smootherName, *laplaceOperator.diagonalValues_, smootherTmp1, smootherTmp2 );
parameters, comm, laplaceOperator, smootherName, *laplaceOperator.getInverseDiagonalValues(), smTmp1, smTmp2, smTmp3 );
WALBERLA_LOG_INFO_ON_ROOT( "created smoother " << smootherName );
auto coarseGridSolver = std::make_shared< CGSolver< P1ElementwiseLaplaceOperator > >(
......
apps/hdgp1/SetupScenario.hpp
View file @ 5bfaa2db
......@@ -42,21 +42,6 @@ namespace setupScenarios {
using CommMode = communication::BufferedCommunicator::LocalCommunicationMode;
template < typename OperatorType, typename FunctionType >
real_t calculateRadius( const OperatorType& op, const FunctionType& diagonal, FunctionType& fct, const uint_t& level )
{
return chebyshev::estimateRadius( op, diagonal, level, 100, fct.getStorage(), fct );
}
template < typename OperatorType >
real_t calculateRadius( const OperatorType& op,
const hdgp1::HDGP1Function< real_t >& diagonal,
hdgp1::HDGP1Function< real_t >& fct,
const uint_t& level )
{
return chebyshev::estimateRadius( op, diagonal, level, 100, fct.getStorage(), fct, fct.getLocalLevelInfo() );
}
template < typename OperatorType, typename FunctionType = typename OperatorType::srcType >
std::shared_ptr< Solver< OperatorType > > setupSpecialSmoother( const walberla::Config::BlockHandle&,
const CommMode&,
......@@ -64,6 +49,7 @@ std::shared_ptr< Solver< OperatorType > > setupSpecialSmoother( const walberla::
std::string&,
FunctionType&,
FunctionType&,
FunctionType&,
FunctionType& )
{
return nullptr;
......@@ -74,18 +60,6 @@ real_t calculateOptimalSmootherWeight( const real_t& radius )
return 1. / ( 0.5 * ( 1.2 * radius + 0.3 * radius ) );
}
template < typename OperatorType >
real_t calculateOptimalSmootherWeight( const OperatorType& op,
hdgp1::HDGP1Function< real_t >& diagonal,
hdgp1::HDGP1Function< real_t >& vec,
const uint_t& level )
{
vec.interpolate( 1., level, All, UpdateMacroEdge );
vec.interpolate( 2., level, All, UpdateMacroFace );
const auto radius = calculateRadius( op, diagonal, vec, level );
return calculateOptimalSmootherWeight( radius );
}
const auto curveConstantInnerValue = []( real_t ) -> real_t { return 0; };
const auto curveConstantOuterValue = []( real_t ) -> real_t { return 1; };
......@@ -104,7 +78,7 @@ enum class OmegaShape
inline bool operator>>( std::istringstream& stream, OmegaShape& shapeEnum )
{
const std::string shape{stream.str()};
const std::string shape{ stream.str() };
if ( shape == "constant_inner" )
shapeEnum = OmegaShape::constant_inner;
else if ( shape == "constant_outer" )
......@@ -142,7 +116,7 @@ enum class EdgeSmootherType
inline bool operator>>( std::istringstream& stream, EdgeSmootherType& edgeEnum )
{
const std::string shape{stream.str()};
const std::string shape{ stream.str() };
if ( shape == "diagonal" )
edgeEnum = EdgeSmootherType::diagonal;
else if ( shape == "block_pairs" )
......@@ -191,14 +165,15 @@ std::shared_ptr< Solver< OperatorType > > setupSpecialSmoother( const walberla::
std::string& smootherName,
hdgp1::HDGP1Function< real_t >& d,
hdgp1::HDGP1Function< real_t >& tmp1,
hdgp1::HDGP1Function< real_t >& tmp2 )
hdgp1::HDGP1Function< real_t >& tmp2,
hdgp1::HDGP1Function< real_t >& tmp3 )
{
const auto minLevel = op.getMinLevel();
const auto maxLevel = op.getMaxLevel();
const auto storage = op.getStorage();
const auto localLevelInfo = d.getLocalLevelInfo();
const std::string smootherType{parameters.getParameter< std::string >( "smootherType" )};
const std::string smootherType{ parameters.getParameter< std::string >( "smootherType" ) };
if ( smootherType == "jacobi_boundary" )
{
......@@ -222,7 +197,7 @@ std::shared_ptr< Solver< OperatorType > > setupSpecialSmoother( const walberla::
const auto curve = toCurve( curveShape );
hdgp1::HDGP1Function< real_t > diagonal( "diagonal", storage, minLevel, maxLevel, localLevelInfo );
op.computeDiagonalOperatorValues( true, diagonal );
op.computeInverseDiagonalOperatorValues( diagonal );
UpdateMacroPrimitive macroPrimitiveFlag = UpdateMacroVertex | UpdateMacroEdge;
if ( !determineOmegaByEdgeRestriction )
......@@ -245,7 +220,7 @@ std::shared_ptr< Solver< OperatorType > > setupSpecialSmoother( const walberla::
}
else
{
dst.multElementwise( {tmp2, diagonal}, level, flag, macroPrimitiveFlag );
dst.multElementwise( { tmp2, diagonal }, level, flag, macroPrimitiveFlag );
dst.interpolate( 0, level, flag, invert( macroPrimitiveFlag ) );
}
} );
......@@ -289,7 +264,7 @@ std::shared_ptr< Solver< OperatorType > > setupSpecialSmoother( const walberla::
const auto getOmega = [&]( const uint_t& level, real_t& radius, real_t& omegaValue ) {
const DoFType flag = DoFType::Inner | DoFType::NeumannBoundary;
utils::fillRandomly( tmp1, level, flag );
radius = estimateSpectralRadiusWithPowerIteration( iterMatrix, tmp1, 100, storage, level, flag, localLevelInfo );
radius = estimateSpectralRadiusWithPowerIteration( iterMatrix, tmp1, tmp3, 100, storage, level, flag );
if ( useLargestOmega )
omegaValue = 2. / radius;
else
......@@ -324,7 +299,7 @@ std::shared_ptr< Solver< OperatorType > > setupSpecialSmoother( const walberla::
WALBERLA_CHECK( pr == faceSmootherFlag );
WALBERLA_CHECK( updateType == Replace, "only replace supported" );
// apply omega and the inverse diagonal on the inner dofs
dst.multElementwise( {faceOmega, diagonal, src}, level, flag, faceSmootherFlag );
dst.multElementwise( { faceOmega, diagonal, src }, level, flag, faceSmootherFlag );
} );
UpdateMacroPrimitive edgeSmootherFlag = UpdateMacroEdge;
......@@ -345,10 +320,10 @@ std::shared_ptr< Solver< OperatorType > > setupSpecialSmoother( const walberla::
}
else
{
dst.multElementwise( {src, diagonal}, level, flag, edgeSmootherFlag );
dst.multElementwise( { src, diagonal }, level, flag, edgeSmootherFlag );
}
// scale with omega
dst.multElementwise( {edgeOmega, dst}, level, flag, edgeSmootherFlag );
dst.multElementwise( { edgeOmega, dst }, level, flag, edgeSmootherFlag );
} );
auto smootherVertices = hdgp1::opalg::wrap(
......@@ -356,7 +331,7 @@ std::shared_ptr< Solver< OperatorType > > setupSpecialSmoother( const walberla::
WALBERLA_CHECK( pr == UpdateMacroVertex );
WALBERLA_CHECK( updateType == Replace, "only replace supported" );
// apply omega and the inverse diagonal on the vertices
dst.multElementwise( {vertexOmega, diagonal, src}, level, flag, UpdateMacroVertex );
dst.multElementwise( { vertexOmega, diagonal, src }, level, flag, UpdateMacroVertex );
} );
auto smoother =
......@@ -405,31 +380,32 @@ std::shared_ptr< Solver< OperatorType > > setupSmoother( const walberla::Config:
std::string& smootherName,
typename OperatorType::srcType& diagonal,
typename OperatorType::srcType& tmp1,
typename OperatorType::srcType& tmp2 )
typename OperatorType::srcType& tmp2,
typename OperatorType::srcType& tmp3 )
{
const auto specialSmoother = setupSpecialSmoother( parameters, commMode, op, smootherName, diagonal, tmp1, tmp2 );
const auto specialSmoother = setupSpecialSmoother( parameters, commMode, op, smootherName, diagonal, tmp1, tmp2, tmp3 );
if ( specialSmoother != nullptr )
return specialSmoother;
const std::string smootherType{parameters.getParameter< std::string >( "smootherType" )};
const std::string smootherType{ parameters.getParameter< std::string >( "smootherType" ) };
// general smoothers
if ( smootherType == "richardson" )
{
auto smoother = std::make_shared< DampedRichardsonSmoother< OperatorType > >( tmp1 );
smoother->setLocalCommunicationMode( commMode );
const auto spectralRadius = calculateRadius( op, diagonal, tmp1, op.getMinLevel() );
const auto spectralRadius = chebyshev::estimateRadius( op, op.getMinLevel(), 100, tmp1.getStorage(), tmp1, tmp2 );
smoother->setWeight( 1. / spectralRadius );
smootherName = smootherType;
return smoother;
}
else if ( smootherType == "chebyshev" )
{
auto smoother = std::make_shared< ChebyshevSmoother< OperatorType > >( diagonal, tmp1, tmp2 );
auto smoother = std::make_shared< ChebyshevSmoother< OperatorType > >( tmp1, tmp2 );
smoother->setLocalCommunicationMode( commMode );
tmp1.interpolate( 0, op.getMinLevel(), DoFType::DirichletBoundary );
const auto spectralRadius = calculateRadius( op, diagonal, tmp1, op.getMinLevel() );
const auto spectralRadius = chebyshev::estimateRadius( op, op.getMinLevel(), 100, tmp1.getStorage(), tmp1, tmp2 );
const uint_t order = parameters.getOneBlock( "Chebyshev" ).getParameter< uint_t >( "order" );
smoother->setupCoefficients( order, spectralRadius );
std::stringstream smootherNameStream{};
......@@ -486,7 +462,7 @@ enum class ProblemType
inline bool operator>>( std::istringstream& stream, ProblemType& problemEnum )
{
const std::string problemString{stream.str()};
const std::string problemString{ stream.str() };
if ( problemString == "sines" )
problemEnum = ProblemType ::sines;
else if ( problemString == "zero_rhs" )
......
apps/hdgp1/Utils.hpp
View file @ 5bfaa2db
......@@ -113,7 +113,7 @@ void setupDiagonal( OperatorType& op )
{
for ( uint_t level = op.getMinLevel(); level <= op.getMaxLevel(); level += 1 )
{
op.computeDiagonalOperatorValues( level, true );
op.computeInverseDiagonalOperatorValues();
}
}
......
src/hyteg/hdgp1functionspace/HDGP1VariableOperator.cpp
View file @ 5bfaa2db
......@@ -100,7 +100,10 @@ HDGP1VariableOperator< FaceForm, MacroEdgeForm, Diagonal, Lumped, InvertDiagonal
const uint_t& maxLevel,
const std::shared_ptr< LocalLevelInfo >& localLevelInfo )
: Operator( storage, minLevel, maxLevel )
, diagonalValues_( new HDGP1Function< real_t >( "inverse diagonal", storage, minLevel, maxLevel, localLevelInfo ) )
// TODO: initialize these vectors when needed.
, diagonalValues_( nullptr )
, inverseDiagonalValues_( nullptr )
, localLevelInfo_( localLevelInfo )
, calculatePenaltyEdgeLength( calculateMinEdgeLength )
{}
......@@ -854,30 +857,42 @@ void HDGP1VariableOperator< FaceForm, MacroEdgeForm, Diagonal, Lumped, InvertDia
}
template < typename FaceForm, typename MacroEdgeForm, bool Diagonal, bool Lumped, bool InvertDiagonal >
void HDGP1VariableOperator< FaceForm, MacroEdgeForm, Diagonal, Lumped, InvertDiagonal >::computeDiagonalOperatorValues(
const bool& invert )
void HDGP1VariableOperator< FaceForm, MacroEdgeForm, Diagonal, Lumped, InvertDiagonal >::computeDiagonalOperatorValues()
{
computeDiagonalOperatorValues( invert, *diagonalValues_ );
if (diagonalValues_ == nullptr)
diagonalValues_ = std::make_shared< HDGP1Function< real_t > >( "diagonal", getStorage(), getMinLevel(), getMaxLevel(), localLevelInfo_ );
computeDiagonalOperatorValues( *diagonalValues_ );
}
template < typename FaceForm, typename MacroEdgeForm, bool Diagonal, bool Lumped, bool InvertDiagonal >
void HDGP1VariableOperator< FaceForm, MacroEdgeForm, Diagonal, Lumped, InvertDiagonal >::computeDiagonalOperatorValues(
const uint_t& level,
const bool& invert )
void HDGP1VariableOperator< FaceForm, MacroEdgeForm, Diagonal, Lumped, InvertDiagonal >::computeInverseDiagonalOperatorValues()
{
computeDiagonalOperatorValues( level, invert, *diagonalValues_ );
if (inverseDiagonalValues_ == nullptr)
inverseDiagonalValues_ = std::make_shared< HDGP1Function< real_t > >( "inverse diagonal", getStorage(), getMinLevel(), getMaxLevel(), localLevelInfo_ );
computeInverseDiagonalOperatorValues( *inverseDiagonalValues_ );
}
template < typename FaceForm, typename MacroEdgeForm, bool Diagonal, bool Lumped, bool InvertDiagonal >
void HDGP1VariableOperator< FaceForm, MacroEdgeForm, Diagonal, Lumped, InvertDiagonal >::computeDiagonalOperatorValues(
const bool& invert,
HDGP1Function< real_t >& dst ) const
void HDGP1VariableOperator< FaceForm, MacroEdgeForm, Diagonal, Lumped, InvertDiagonal >::computeDiagonalOperatorValues( HDGP1Function< real_t >& dst ) const
{
WALBERLA_ASSERT( !InvertDiagonal, "InvertDiagonal not implemented yet." );
for ( uint_t levelToFill = minLevel_; levelToFill <= maxLevel_; levelToFill += 1 )
{
computeDiagonalOperatorValues( levelToFill, false, dst );
}
}
template < typename FaceForm, typename MacroEdgeForm, bool Diagonal, bool Lumped, bool InvertDiagonal >
void HDGP1VariableOperator< FaceForm, MacroEdgeForm, Diagonal, Lumped, InvertDiagonal >::computeInverseDiagonalOperatorValues( HDGP1Function< real_t >& dst ) const
{
WALBERLA_ASSERT( !InvertDiagonal, "InvertDiagonal not implemented yet." );
for ( uint_t levelToFill = minLevel_; levelToFill <= maxLevel_; levelToFill += 1 )
{
computeDiagonalOperatorValues( levelToFill, invert, dst );
computeDiagonalOperatorValues( levelToFill, true, dst );
}
}
......@@ -1268,7 +1283,7 @@ void HDGP1VariableOperator< FaceForm, MacroEdgeForm, Diagonal, Lumped, InvertDia
}
}
const real_t* const diagonalData = face.getData( diagonalValues_->getFaceDataID() )->getPointer( level );
const real_t* const diagonalData = face.getData( getInverseDiagonalValues()->getFaceDataID() )->getPointer( level );
const std::array< uint_t, 3 > vIndices{
0, levelinfo::num_microvertices_per_edge( localLevel ) - 1, levelinfo::num_microvertices_per_face( localLevel ) - 1 };
......@@ -1457,7 +1472,7 @@ void HDGP1VariableOperator< FaceForm, MacroEdgeForm, Diagonal, Lumped, InvertDia
}
}
const real_t* const diagonalData = face.getData( diagonalValues_->getFaceDataID() )->getPointer( level );
const real_t* const diagonalData = face.getData( getInverseDiagonalValues()->getFaceDataID() )->getPointer( level );
const std::array< uint_t, 3 > vIndices{
0, levelinfo::num_microvertices_per_edge( localLevel ) - 1, levelinfo::num_microvertices_per_face( localLevel ) - 1 };
......
src/hyteg/hdgp1functionspace/HDGP1VariableOperator.hpp
View file @ 5bfaa2db
......@@ -63,12 +63,26 @@ class HDGP1VariableOperator : public Operator< HDGP1Function< real_t >, HDGP1Fun
UpdateType updateType = Replace,
UpdateMacroPrimitive primitive = UpdateAllPrimitives ) const;
std::shared_ptr< HDGP1Function< real_t > > diagonalValues_;
void computeDiagonalOperatorValues( const bool& invert );
void computeDiagonalOperatorValues( const uint_t& level, const bool& invert );
void computeDiagonalOperatorValues( const bool& invert, HDGP1Function< real_t >& dst ) const;
void computeDiagonalOperatorValues( const uint_t& level, const bool& invert, HDGP1Function< real_t >& dst ) const;
void computeDiagonalOperatorValues();
void computeInverseDiagonalOperatorValues();
void computeDiagonalOperatorValues( HDGP1Function< real_t >& dst ) const;
void computeInverseDiagonalOperatorValues( HDGP1Function< real_t >& dst ) const;
std::shared_ptr< HDGP1Function< real_t > > getDiagonalValues() const
{
WALBERLA_CHECK_NOT_NULLPTR(
diagonalValues_,
"Diagonal values have not been assembled, call computeDiagonalOperatorValues() to set up this function." )
return diagonalValues_;
};
std::shared_ptr< HDGP1Function< real_t > > getInverseDiagonalValues() const
{
WALBERLA_CHECK_NOT_NULLPTR(
inverseDiagonalValues_,
"Inverse diagonal values have not been assembled, call computeInverseDiagonalOperatorValues() to set up this function." )
return inverseDiagonalValues_;
};
void applyInverseBlockDiagonal( const HDGP1Function< real_t >& src,
const HDGP1Function< real_t >& dst,
......@@ -89,6 +103,8 @@ class HDGP1VariableOperator : public Operator< HDGP1Function< real_t >, HDGP1Fun
void useMinTangentialEdgeLengthAsPenaltyEdgeLength();
private:
void computeDiagonalOperatorValues( const uint_t& level, const bool& invert, HDGP1Function< real_t >& dst ) const;
template < bool Transpose >
void applyGeneric( const HDGP1Function< real_t >& src,
const HDGP1Function< real_t >& dst,
......@@ -136,6 +152,11 @@ class HDGP1VariableOperator : public Operator< HDGP1Function< real_t >, HDGP1Fun
FaceForm faceForm;
MacroEdgeForm edgeForm;
std::shared_ptr< HDGP1Function< real_t > > diagonalValues_;
std::shared_ptr< HDGP1Function< real_t > > inverseDiagonalValues_;
const std::shared_ptr< LocalLevelInfo >& localLevelInfo_;
std::function< real_t( std::array< Point3D, 4 > ) > calculatePenaltyEdgeLength;
};
......
src/hyteg/numerictools/SpectrumEstimation.hpp
View file @ 5bfaa2db
......@@ -37,16 +37,6 @@ namespace hyteg {
// =================================================================================================
template < typename OperatorType >
real_t estimateSpectralRadiusWithPowerIteration( OperatorType& op,
typename OperatorType::srcType& itrVec,
const uint_t numIts,
const std::shared_ptr< PrimitiveStorage >& storage,
const uint_t level )
{
return estimateSpectralRadiusWithPowerIteration( op, itrVec, numIts, storage, level, DoFType::All );
}
/// \brief Compute an estimate for the spectral radius of an operator using the Power Method.
///
/// The estimate for the spetral radius is computed using a simple Power iteration.
......@@ -69,7 +59,7 @@ namespace hyteg {
const uint_t numIts,
const std::shared_ptr<PrimitiveStorage> &storage,
const uint_t level,
const DoFType flag ) {
const DoFType flag = DoFType::All ) {
// scale input vector
real_t norm = std::sqrt( itrVec.dotGlobal( itrVec, level, flag ) );
......
src/hyteg/solvers/ChebyshevSmoother.hpp
View file @ 5bfaa2db
......@@ -41,6 +41,14 @@ class ChebyshevSmoother : public Solver< OperatorType >
, flag_( Inner | NeumannBoundary )
{}
ChebyshevSmoother( const FunctionType& tmp1, const FunctionType& tmp2 )
: coefficients{}
, tmp1_( tmp1 )
, tmp2_( tmp2 )
, flag_( Inner | NeumannBoundary )
{}
void setLocalCommunicationMode( const communication::BufferedCommunicator::LocalCommunicationMode& localCommunicationMode )
{