Commit 513aead8 authored by Nils Kohl's avatar Nils Kohl 🌝
Browse files

Fixing tests and project normal operators.

parent fbff74ae
......@@ -22,6 +22,7 @@
#include "hyteg/composites/P1StokesBlockPreconditioner.hpp"
#include "hyteg/composites/P1StokesFunction.hpp"
#include "hyteg/composites/StokesOperatorTraits.hpp"
#include "hyteg/operators/VectorLaplaceOperator.hpp"
#include "hyteg/p1functionspace/P1ConstantOperator.hpp"
#include "hyteg/p1functionspace/P1ScalarToP1VectorOperator.hpp"
#include "hyteg/p1functionspace/P1VectorToP1ScalarOperator.hpp"
......@@ -31,9 +32,10 @@ namespace hyteg {
class P1StokesOperator : public Operator< P1StokesFunction< real_t >, P1StokesFunction< real_t > >
{
public:
typedef P1ConstantLaplaceOperator VelocityOperator_T;
typedef P1ConstantLaplaceOperator PressureOperator_T;
typedef P1StokesBlockPreconditioner BlockPreconditioner_T;
typedef P1ConstantVectorLaplaceOperator VelocityBlockOperator_T;
typedef P1ConstantLaplaceOperator VelocityOperator_T;
typedef P1ConstantLaplaceOperator PressureOperator_T;
typedef P1StokesBlockPreconditioner BlockPreconditioner_T;
P1StokesOperator( const std::shared_ptr< PrimitiveStorage >& storage, size_t minLevel, size_t maxLevel )
: Operator( storage, minLevel, maxLevel )
......
......@@ -33,13 +33,21 @@ namespace hyteg {
class P2P1TaylorHoodStokesOperator : public Operator< P2P1TaylorHoodFunction< real_t >, P2P1TaylorHoodFunction< real_t > >
{
public:
typedef P2ConstantVectorLaplaceOperator VelocityBlockOperator_T;
typedef P2ConstantLaplaceOperator VelocityOperator_T;
typedef P2P1TaylorHoodStokesBlockPreconditioner BlockPreconditioner_T;
P2P1TaylorHoodStokesOperator( const std::shared_ptr< PrimitiveStorage >& storage, size_t minLevel, size_t maxLevel )
: Operator( storage, minLevel, maxLevel )
, A( storage, minLevel, maxLevel )
, Lapl( storage, minLevel, maxLevel )
, div_x( storage, minLevel, maxLevel )
, div_y( storage, minLevel, maxLevel )
, div_z( storage, minLevel, maxLevel )
, div( storage, minLevel, maxLevel )
, divT_x( storage, minLevel, maxLevel )
, divT_y( storage, minLevel, maxLevel )
, divT_z( storage, minLevel, maxLevel )
, divT( storage, minLevel, maxLevel )
, pspg_( storage, minLevel, maxLevel )
, pspg_inv_diag_( storage, minLevel, maxLevel )
......@@ -67,9 +75,18 @@ class P2P1TaylorHoodStokesOperator : public Operator< P2P1TaylorHoodFunction< re
div.toMatrix( mat, src.uvw, dst.p, level, flag );
}
P2ConstantVectorLaplaceOperator Lapl;
P2ToP1ConstantDivOperator div;
P1ToP2ConstantDivTOperator divT;
VelocityBlockOperator_T Lapl;
P2ToP1ConstantDivOperator div;
P1ToP2ConstantDivTOperator divT;
// currently, need these for being able to integrate into the framework until the block operator is completed
P2ConstantLaplaceOperator A;
P2ToP1ConstantDivxOperator div_x;
P2ToP1ConstantDivyOperator div_y;
P2ToP1ConstantDivzOperator div_z;
P1ToP2ConstantDivTxOperator divT_x;
P1ToP2ConstantDivTyOperator divT_y;
P1ToP2ConstantDivTzOperator divT_z;
// this operator is needed in the uzawa smoother
P1PSPGOperator pspg_;
......
......@@ -101,8 +101,8 @@ class StrongFreeSlipWrapper : public Operator< typename OpType::srcType, typenam
/// \param flag determines on which primitives this operator is assembled
///
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const typename OpType::srcType::template FunctionType< PetscInt >& numeratorSrc,
const typename OpType::dstType::template FunctionType< PetscInt >& numeratorDst,
const typename OpType::srcType::template FunctionType< matIdx_t >& numeratorSrc,
const typename OpType::dstType::template FunctionType< matIdx_t >& numeratorDst,
uint_t level,
DoFType flag ) const
{
......
......@@ -20,6 +20,7 @@
#include <typeinfo>
#include "hyteg/elementwiseoperators/P2P1ElementwiseBlendingStokesOperator.hpp"
#include "hyteg/functions/FunctionProperties.hpp"
#include "hyteg/indexing/CouplingCount.hpp"
#include "hyteg/p1functionspace/P1Function.hpp"
......@@ -46,8 +47,8 @@ uint_t getNumberOfGlobalDoFCouplings( const opType& oper, uint_t level )
{
uint_t nCouplings = 0;
typedef Operator< P1Function< real_t >, P1Function< real_t > > P1ScalarOp;
typedef Operator< P2Function< real_t >, P2Function< real_t > > P2ScalarOp;
typedef Operator< P1Function< real_t >, P1Function< real_t > > P1ScalarOp;
typedef Operator< P2Function< real_t >, P2Function< real_t > > P2ScalarOp;
typedef Operator< EdgeDoFFunction< real_t >, EdgeDoFFunction< real_t > > EdgeDoFScalarOp;
// Scalar P1 operators
......
......@@ -20,11 +20,11 @@
#include "P1ProjectNormalOperator.hpp"
#include "hyteg/communication/Syncing.hpp"
#include "hyteg/p1functionspace/freeslip/VertexDoFProjectNormal.hpp"
#include "hyteg/p1functionspace/VertexDoFMacroVertex.hpp"
#include "hyteg/p1functionspace/VertexDoFMacroCell.hpp"
#include "hyteg/p1functionspace/VertexDoFMacroEdge.hpp"
#include "hyteg/p1functionspace/VertexDoFMacroFace.hpp"
#include "hyteg/p1functionspace/VertexDoFMacroCell.hpp"
#include "hyteg/p1functionspace/VertexDoFMacroVertex.hpp"
#include "hyteg/p1functionspace/freeslip/VertexDoFProjectNormal.hpp"
namespace hyteg {
......@@ -266,14 +266,12 @@ void P1ProjectNormalOperator::project( const P1StokesFunction< real_t >& dst, si
this->stopTiming( "Project" );
}
#ifdef HYTEG_BUILD_WITH_PETSC
void P1ProjectNormalOperator::assembleLocalMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P1Function< PetscInt >& numU,
const P1Function< PetscInt >& numV,
const P1Function< PetscInt >& numW,
uint_t level,
DoFType flag ) const
void P1ProjectNormalOperator::toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P1Function< matIdx_t >& numU,
const P1Function< matIdx_t >& numV,
const P1Function< matIdx_t >& numW,
uint_t level,
DoFType flag ) const
{
communication::syncFunctionBetweenPrimitives( numU, level );
communication::syncFunctionBetweenPrimitives( numV, level );
......@@ -373,12 +371,7 @@ void P1ProjectNormalOperator::assembleLocalMatrix( const std::shared_ptr< Sparse
}
}
}
}
}
#endif
} // namespace hyteg
......@@ -20,16 +20,16 @@
#pragma once
#include "hyteg/HytegDefinitions.hpp"
#include "hyteg/sparseassembly/SparseMatrixProxy.hpp"
#include "hyteg/operators/Operator.hpp"
#include "hyteg/composites//P1StokesFunction.hpp"
#include "hyteg/operators/Operator.hpp"
#include "hyteg/petsc/PETScWrapper.hpp"
#include "hyteg/sparseassembly/SparseMatrixProxy.hpp"
namespace hyteg {
using walberla::real_t;
class P1ProjectNormalOperator : public Operator< P1Function< real_t >, P1Function< real_t > >
class P1ProjectNormalOperator : public Operator< P1VectorFunction< real_t >, P1VectorFunction< real_t > >
{
public:
P1ProjectNormalOperator( const std::shared_ptr< PrimitiveStorage >& storage,
......@@ -45,9 +45,17 @@ class P1ProjectNormalOperator : public Operator< P1Function< real_t >, P1Functio
size_t level,
DoFType flag ) const;
void project( const P1VectorFunction< real_t >& dst, size_t level, DoFType flag ) const;
void project( const P1StokesFunction< real_t >& dst, size_t level, DoFType flag ) const;
#ifdef HYTEG_BUILD_WITH_PETSC
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P1Function< matIdx_t >& numU,
const P1Function< matIdx_t >& numV,
const P1Function< matIdx_t >& numW,
uint_t level,
DoFType flag ) const;
/// Assemble operator as sparse matrix
///
/// \param mat a sparse matrix proxy
......@@ -57,13 +65,15 @@ class P1ProjectNormalOperator : public Operator< P1Function< real_t >, P1Functio
/// \param level level in mesh hierarchy for which local operator is to be assembled
/// \param flag determines on which primitives this operator is assembled
///
void assembleLocalMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P1Function< PetscInt >& numU,
const P1Function< PetscInt >& numV,
const P1Function< PetscInt >& numW,
uint_t level,
DoFType flag ) const;
#endif
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P1VectorFunction< matIdx_t >& src,
const P1VectorFunction< matIdx_t >& dst,
uint_t level,
DoFType flag ) const override
{
WALBERLA_UNUSED( dst );
toMatrix( mat, src[0], src[1], src[2], level, flag );
}
private:
const std::function< void( const Point3D&, Point3D& ) > normal_function_;
......
......@@ -51,31 +51,26 @@ void P2ProjectNormalOperator::project( const P2P1TaylorHoodFunction< real_t >& d
project( dst.uvw[0], dst.uvw[1], dst.uvw[idx], level, flag );
}
#ifdef HYTEG_BUILD_WITH_PETSC
void P2ProjectNormalOperator::assembleLocalMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P2Function< PetscInt >& numU,
const P2Function< PetscInt >& numV,
const P2Function< PetscInt >& numW,
uint_t level,
DoFType flag )
void P2ProjectNormalOperator::toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P2Function< matIdx_t >& numU,
const P2Function< matIdx_t >& numV,
const P2Function< matIdx_t >& numW,
uint_t level,
DoFType flag ) const
{
p1Operator.assembleLocalMatrix(
mat, numU.getVertexDoFFunction(), numV.getVertexDoFFunction(), numW.getVertexDoFFunction(), level, flag );
p1Operator.toMatrix( mat, numU.getVertexDoFFunction(), numV.getVertexDoFFunction(), numW.getVertexDoFFunction(), level, flag );
edgeDoFOperator.assembleLocalMatrix(
mat, numU.getEdgeDoFFunction(), numV.getEdgeDoFFunction(), numW.getEdgeDoFFunction(), level, flag );
}
void P2ProjectNormalOperator::assembleLocalMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P2VectorFunction< PetscInt >& num,
uint_t level,
DoFType flag )
void P2ProjectNormalOperator::toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P2VectorFunction< matIdx_t >& num,
uint_t level,
DoFType flag ) const
{
// UGLY FIX (for 2D the 3rd component function is not accessed later on anyway!)
uint_t idx = num.getDimension() == 2 ? 0 : 2;
assembleLocalMatrix( mat, num[0], num[1], num[idx], level, flag );
toMatrix( mat, num[0], num[1], num[idx], level, flag );
}
#endif
} // namespace hyteg
......@@ -19,11 +19,11 @@
*/
#pragma once
#include "hyteg/operators/Operator.hpp"
#include "hyteg/composites/P2P1TaylorHoodFunction.hpp"
#include "hyteg/composites/P2P1TaylorHoodStokesOperator.hpp"
#include "hyteg/composites/petsc/P2P1TaylorHoodPetsc.hpp"
#include "hyteg/edgedofspace/EdgeDoFProjectNormalOperator.hpp"
#include "hyteg/operators/Operator.hpp"
#include "hyteg/p1functionspace/P1ProjectNormalOperator.hpp"
#include "hyteg/p2functionspace/P2Function.hpp"
......@@ -31,7 +31,7 @@ namespace hyteg {
using walberla::real_t;
class P2ProjectNormalOperator : public Operator< P2Function< real_t >, P2Function< real_t > >
class P2ProjectNormalOperator : public Operator< P2VectorFunction< real_t >, P2VectorFunction< real_t > >
{
public:
P2ProjectNormalOperator( const std::shared_ptr< PrimitiveStorage >& storage,
......@@ -47,9 +47,13 @@ class P2ProjectNormalOperator : public Operator< P2Function< real_t >, P2Functio
size_t level,
DoFType flag ) const;
void project( const P2VectorFunction< real_t >& dst, size_t level, DoFType flag ) const
{
project( dst[0], dst[1], dst[2], level, flag );
}
void project( const P2P1TaylorHoodFunction< real_t >& dst, size_t level, DoFType flag ) const;
#ifdef HYTEG_BUILD_WITH_PETSC
/// Assemble operator as sparse matrix
///
/// \param mat a sparse matrix proxy
......@@ -59,12 +63,17 @@ class P2ProjectNormalOperator : public Operator< P2Function< real_t >, P2Functio
/// \param level level in mesh hierarchy for which local operator is to be assembled
/// \param flag determines on which primitives this operator is assembled
///
void assembleLocalMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P2Function< PetscInt >& numU,
const P2Function< PetscInt >& numV,
const P2Function< PetscInt >& numW,
uint_t level,
DoFType flag );
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P2Function< matIdx_t >& numU,
const P2Function< matIdx_t >& numV,
const P2Function< matIdx_t >& numW,
uint_t level,
DoFType flag ) const;
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P2VectorFunction< matIdx_t >& num,
uint_t level,
DoFType flag ) const;
/// Assemble operator as sparse matrix
///
......@@ -73,16 +82,19 @@ class P2ProjectNormalOperator : public Operator< P2Function< real_t >, P2Functio
/// \param level level in mesh hierarchy for which local operator is to be assembled
/// \param flag determines on which primitives this operator is assembled
///
void assembleLocalMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P2VectorFunction< PetscInt >& num,
uint_t level,
DoFType flag );
#endif
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P2VectorFunction< matIdx_t >& src,
const P2VectorFunction< matIdx_t >& dst,
uint_t level,
DoFType flag ) const override
{
WALBERLA_UNUSED( dst );
toMatrix( mat, src, level, flag );
}
private:
P1ProjectNormalOperator p1Operator;
EdgeDoFProjectNormalOperator edgeDoFOperator;
};
} // namespace hyteg
......@@ -200,8 +200,7 @@ class PETScBlockPreconditionedStokesSolver : public Solver< OperatorType >
case 1:
PCSetType( pc_u, PCJACOBI );
break;
case 3:
{
case 3: {
PCSetType( pc_u, PCHYPRE );
break;
}
......@@ -225,9 +224,7 @@ class PETScBlockPreconditionedStokesSolver : public Solver< OperatorType >
default:
WALBERLA_ABORT( "Invalid pressure preconditioner for PETSc block prec MinRes solver." );
break;
}
}
timer.end();
......@@ -302,7 +299,8 @@ class PETScBlockPreconditionedStokesSolver : public Solver< OperatorType >
uint_t level,
const P2P1TaylorHoodFunction< PetscInt >& numerator )
{
for ( auto dof : FunctionIterator< vertexdof::VertexDoFFunction< PetscInt > >( numerator.uvw[0].getVertexDoFFunction(), level ) )
for ( auto dof :
FunctionIterator< vertexdof::VertexDoFFunction< PetscInt > >( numerator.uvw[0].getVertexDoFFunction(), level ) )
{
velocityIndices.push_back( dof.value() );
}
......@@ -310,7 +308,8 @@ class PETScBlockPreconditionedStokesSolver : public Solver< OperatorType >
{
velocityIndices.push_back( dof.value() );
}
for ( auto dof : FunctionIterator< vertexdof::VertexDoFFunction< PetscInt > >( numerator.uvw[1].getVertexDoFFunction(), level ) )
for ( auto dof :
FunctionIterator< vertexdof::VertexDoFFunction< PetscInt > >( numerator.uvw[1].getVertexDoFFunction(), level ) )
{
velocityIndices.push_back( dof.value() );
}
......@@ -336,15 +335,15 @@ class PETScBlockPreconditionedStokesSolver : public Solver< OperatorType >
}
}
uint_t allocatedLevel_;
MPI_Comm petscCommunicator_;
typename OperatorType::srcType::template FunctionType< PetscInt > num;
PETScSparseMatrix< OperatorType, OperatorType::srcType::template FunctionType > Amat;
PETScSparseMatrix< OperatorType, OperatorType::srcType::template FunctionType > AmatNonEliminatedBC;
PETScSparseMatrix< BlockPreconditioner_T, BlockPreconditioner_T::srcType::template FunctionType > Pmat;
PETScVector< typename FunctionType::valueType, OperatorType::srcType::template FunctionType > xVec;
PETScVector< typename FunctionType::valueType, OperatorType::srcType::template FunctionType > bVec;
PETScVector< typename FunctionType::valueType, OperatorType::srcType::template FunctionType > nullspaceVec_;
uint_t allocatedLevel_;
MPI_Comm petscCommunicator_;
typename OperatorType::srcType::template FunctionType< PetscInt > num;
PETScSparseMatrix< OperatorType > Amat;
PETScSparseMatrix< OperatorType > AmatNonEliminatedBC;
PETScSparseMatrix< BlockPreconditioner_T > Pmat;
PETScVector< typename FunctionType::valueType, OperatorType::srcType::template FunctionType > xVec;
PETScVector< typename FunctionType::valueType, OperatorType::srcType::template FunctionType > bVec;
PETScVector< typename FunctionType::valueType, OperatorType::srcType::template FunctionType > nullspaceVec_;
std::shared_ptr< PrimitiveStorage > storage_;
......
......@@ -92,8 +92,8 @@ void exportLinearSystem( OperatorType op,
{
WALBERLA_LOG_INFO_ON_ROOT( " * Converting Operator to PETSc matrix" );
}
PETScSparseMatrix< OperatorType, FunctionType > petscMatrix( localDoFs, globalDoFs, nameMatrix.c_str() );
FunctionType< PetscInt > numerator( "numerator", storage, level, level );
PETScSparseMatrix< OperatorType > petscMatrix( localDoFs, globalDoFs, nameMatrix.c_str() );
FunctionType< PetscInt > numerator( "numerator", storage, level, level );
numerator.enumerate( level );
petscMatrix.createMatrixFromOperator( op, level, numerator );
......
......@@ -80,7 +80,7 @@ void compareCounts( std::shared_ptr< PrimitiveStorage > storage, std::string tag
uint_t nnzHyTeG = indexing::getNumberOfGlobalDoFCouplings( oper, level );
// assemble matrix
PETScSparseMatrix< opType, fType > petscMat( localDoFs, globalDoFs );
PETScSparseMatrix< opType > petscMat( localDoFs, globalDoFs );
petscMat.createMatrixFromOperator( oper, level, enumerator, All );
uint_t nnzPETSc = petscMat.getInfo().getNNZ();
......
......@@ -42,95 +42,97 @@ using walberla::real_t;
namespace hyteg {
void p1PetscApplyTest( const uint_t & level, const std::string & meshFile, const DoFType & location, const real_t & eps )
void p1PetscApplyTest( const uint_t& level, const std::string& meshFile, const DoFType& location, const real_t& eps )
{
WALBERLA_LOG_INFO_ON_ROOT( "level: " << level << ", mesh file: " << meshFile );
WALBERLA_LOG_INFO_ON_ROOT( "level: " << level << ", mesh file: " << meshFile );
MeshInfo meshInfo = hyteg::MeshInfo::fromGmshFile( meshFile );
SetupPrimitiveStorage setupStorage(meshInfo, walberla::uint_c(walberla::mpi::MPIManager::instance()->numProcesses()));
setupStorage.setMeshBoundaryFlagsOnBoundary( 1, 0, true );
loadbalancing::roundRobin( setupStorage );
std::shared_ptr< hyteg::PrimitiveStorage> storage = std::make_shared< hyteg::PrimitiveStorage>(setupStorage);
MeshInfo meshInfo = hyteg::MeshInfo::fromGmshFile( meshFile );
SetupPrimitiveStorage setupStorage( meshInfo, walberla::uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) );
setupStorage.setMeshBoundaryFlagsOnBoundary( 1, 0, true );
loadbalancing::roundRobin( setupStorage );
std::shared_ptr< hyteg::PrimitiveStorage > storage = std::make_shared< hyteg::PrimitiveStorage >( setupStorage );
writeDomainPartitioningVTK( storage, "../../output", "P1PetscApplyTestDomain" );
writeDomainPartitioningVTK( storage, "../../output", "P1PetscApplyTestDomain" );
P1Function< real_t > src ( "src", storage, level, level );
P1Function< real_t > hhgDst ( "hhgDst", storage, level, level );
P1Function< real_t > petscDst ( "petscDst", storage, level, level );
P1Function< real_t > err ( "error", storage, level, level );
P1Function< real_t > ones ( "ones", storage, level, level );
P1Function< PetscInt > numerator( "numerator", storage, level, level );
P1Function< real_t > src( "src", storage, level, level );
P1Function< real_t > hhgDst( "hhgDst", storage, level, level );
P1Function< real_t > petscDst( "petscDst", storage, level, level );
P1Function< real_t > err( "error", storage, level, level );
P1Function< real_t > ones( "ones", storage, level, level );
P1Function< PetscInt > numerator( "numerator", storage, level, level );
std::function<real_t(const hyteg::Point3D&)> zero = [](const hyteg::Point3D&) { return 0.0; };
std::function<real_t(const hyteg::Point3D&)> one = [](const hyteg::Point3D&) { return 1.0; };
std::function<real_t(const hyteg::Point3D&)> rand = []( const hyteg::Point3D & ) { return walberla::math::realRandom<real_t>(); };
std::function<real_t(const hyteg::Point3D&)> srcFunction = []( const hyteg::Point3D & x ) { return x[0] * x[0] * x[0] * x[0] * std::sinh( x[1] ) * std::cos( x[2] ); };
std::function< real_t( const hyteg::Point3D& ) > zero = []( const hyteg::Point3D& ) { return 0.0; };
std::function< real_t( const hyteg::Point3D& ) > one = []( const hyteg::Point3D& ) { return 1.0; };
std::function< real_t( const hyteg::Point3D& ) > rand = []( const hyteg::Point3D& ) {
return walberla::math::realRandom< real_t >();
};
std::function< real_t( const hyteg::Point3D& ) > srcFunction = []( const hyteg::Point3D& x ) {
return x[0] * x[0] * x[0] * x[0] * std::sinh( x[1] ) * std::cos( x[2] );
};
src.interpolate( srcFunction, level, hyteg::All );
hhgDst.interpolate( rand, level, location );
petscDst.interpolate( rand, level, location );
ones.interpolate( one, level, location );
src.interpolate( srcFunction, level, hyteg::All );
hhgDst.interpolate( rand, level, location );
petscDst.interpolate( rand, level, location );
ones.interpolate( one, level, location );
P1ConstantLaplaceOperator L( storage, level, level );
P1ConstantLaplaceOperator L( storage, level, level );
numerator.enumerate( level );
numerator.enumerate( level );
const uint_t globalDoFs = hyteg::numberOfGlobalDoFs< hyteg::P1FunctionTag >( *storage, level );
const uint_t localDoFs = hyteg::numberOfLocalDoFs< hyteg::P1FunctionTag >( *storage, level );
const uint_t globalDoFs = hyteg::numberOfGlobalDoFs< hyteg::P1FunctionTag >( *storage, level );
const uint_t localDoFs = hyteg::numberOfLocalDoFs< hyteg::P1FunctionTag >( *storage, level );
WALBERLA_LOG_INFO_ON_ROOT( "Global DoFs: " << globalDoFs );
WALBERLA_LOG_INFO_ON_ROOT( "Global DoFs: " << globalDoFs );
// HyTeG apply
L.apply( src, hhgDst, level, location );
// HyTeG apply
L.apply( src, hhgDst, level, location );
// PETSc apply
PETScVector< real_t, P1Function > srcPetscVec( localDoFs );
PETScVector< real_t, P1Function > dstPetscVec( localDoFs );
PETScSparseMatrix< P1ConstantLaplaceOperator, P1Function > petscMatrix( localDoFs, globalDoFs );
// PETSc apply
PETScVector< real_t, P1Function > srcPetscVec( localDoFs );
PETScVector< real_t, P1Function > dstPetscVec( localDoFs );
PETScSparseMatrix< P1ConstantLaplaceOperator > petscMatrix( localDoFs, globalDoFs );
srcPetscVec.createVectorFromFunction( src, numerator, level, All );
dstPetscVec.createVectorFromFunction( petscDst, numerator, level, All );
petscMatrix.createMatrixFromOperator( L, level, numerator, All );
srcPetscVec.createVectorFromFunction( src, numerator, level, All );
dstPetscVec.createVectorFromFunction( petscDst, numerator, level, All );
petscMatrix.createMatrixFromOperator( L, level, numerator, All );
WALBERLA_CHECK( petscMatrix.isSymmetric() );
WALBERLA_CHECK( petscMatrix.isSymmetric() );
MatMult( petscMatrix.get(), srcPetscVec.get(), dstPetscVec.get() );
MatMult( petscMatrix.get(), srcPetscVec.get(), dstPetscVec.get() );
dstPetscVec.createFunctionFromVector( petscDst, numerator, level, location );
dstPetscVec.createFunctionFromVector( petscDst, numerator, level, location );
// compare
err.assign( {1.0, -1.0}, {hhgDst, petscDst}, level, location );
const auto absScalarProd = std::abs( err.dotGlobal( ones, level, location ) );
// compare
err.assign( { 1.0, -1.0 }, { hhgDst, petscDst }, level, location );
const auto absScalarProd = std::abs( err.dotGlobal( ones, level, location ) );
WALBERLA_LOG_INFO_ON_ROOT( "Error sum = " << absScalarProd );
WALBERLA_LOG_INFO_ON_ROOT( "Error sum = " << absScalarProd );
// VTK
// VTKOutput vtkOutput( "../../output", "P1PetscApplyTest", storage );
// vtkOutput.add( src );
// vtkOutput.add( hhgDst );
// vtkOutput.add( petscDst );
// vtkOutput.add( err );
// vtkOutput.write( level, 0 );
WALBERLA_CHECK_LESS( absScalarProd, eps );
// VTK
// VTKOutput vtkOutput( "../../output", "P1PetscApplyTest", storage );
// vtkOutput.add( src );
// vtkOutput.add( hhgDst );
// vtkOutput.add( petscDst );
// vtkOutput.add( err );
// vtkOutput.write( level, 0 );
WALBERLA_CHECK_LESS( absScalarProd, eps );
}
}
} // namespace hyteg
int main(int argc, char* argv[])
int main( int argc, char* argv[] )
{
walberla::MPIManager::instance()->initializeMPI( &argc, &argv );
walberla::MPIManager::instance()->useWorldComm();
hyteg::PETScManager petscManager( &argc, &argv );
hyteg::p1PetscApplyTest( 3, "../../data/meshes/quad_4el.msh", hyteg::All, 1.0e-15 );
hyteg::p1PetscApplyTest( 3, "../../data/meshes/annulus_coarse.msh", hyteg::All, 1.7e-13 );
hyteg::p1PetscApplyTest( 3, "../../data/meshes/3D/tet_1el.msh", hyteg::Inner, 7.6e-18 );
hyteg::p1PetscApplyTest( 3, "../../data/meshes/3D/pyramid_2el.msh", hyteg::Inner, 2.0e-16 );
hyteg::p1PetscApplyTest( 3, "../../data/meshes/3D/pyramid_4el.msh", hyteg::Inner, 1.0e-15 );
hyteg::p1PetscApplyTest( 3, "../../data/meshes/3D/regular_octahedron_8el.msh", hyteg::Inner, 1.0e-15 );