Commit a3b8ed20 authored by Marcel Koch's avatar Marcel Koch
Browse files

fixup! ginkgo block matrix test

parent 03be1737
......@@ -17,9 +17,7 @@
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <ginkgo/core/distributed/communicator.hpp>
#include <numeric>
#include "core/DataTypes.h"
#include "core/math/Random.h"
......@@ -28,14 +26,14 @@
#include "hyteg/communication/Syncing.hpp"
#include "hyteg/composites/P2P1TaylorHoodFunction.hpp"
#include "hyteg/composites/P2P1TaylorHoodStokesOperator.hpp"
#include "hyteg/composites/petsc/P2P1TaylorHoodPetsc.hpp"
#include "hyteg/dataexport/VTKOutput.hpp"
#include "hyteg/functions/FunctionIterator.hpp"
#include "hyteg/functions/FunctionTraits.hpp"
#include "hyteg/ginkgo/GinkgoSparseMatrixProxy.hpp"
#include "hyteg/ginkgo/GinkgoUtilities.hpp"
#include "hyteg/ginkgo/GinkgoVectorProxy.hpp"
#include "hyteg/mesh/MeshInfo.hpp"
#include "hyteg/petsc/PETScManager.hpp"
#include "hyteg/petsc/PETScVector.hpp"
#include "hyteg/primitivestorage/PrimitiveStorage.hpp"
#include "hyteg/primitivestorage/SetupPrimitiveStorage.hpp"
#include "hyteg/primitivestorage/Visualization.hpp"
......@@ -45,78 +43,82 @@ using walberla::real_t;
namespace hyteg {
void gatherIndices( std::vector< PetscInt >& velocityIndices,
std::vector< PetscInt >& pressureIndices,
const PrimitiveStorage& storage,
uint_t level,
const P1StokesFunction< PetscInt >& numerator )
void gatherIndices( std::vector< int32_t >& velocityIndices,
std::vector< int32_t >& pressureIndices,
const PrimitiveStorage& storage,
uint_t level,
const P1StokesFunction< int32_t >& numerator )
{
for ( auto dof : FunctionIterator< vertexdof::VertexDoFFunction< PetscInt > >( numerator.uvw[0], level ) )
for ( auto dof : FunctionIterator< vertexdof::VertexDoFFunction< int32_t > >( numerator.uvw[0], level ) )
{
velocityIndices.push_back( dof.value() );
}
for ( auto dof : FunctionIterator< vertexdof::VertexDoFFunction< PetscInt > >( numerator.uvw[1], level ) )
for ( auto dof : FunctionIterator< vertexdof::VertexDoFFunction< int32_t > >( numerator.uvw[1], level ) )
{
velocityIndices.push_back( dof.value() );
}
if ( storage.hasGlobalCells() )
{
for ( auto dof : FunctionIterator< vertexdof::VertexDoFFunction< PetscInt > >( numerator.uvw[2], level ) )
for ( auto dof : FunctionIterator< vertexdof::VertexDoFFunction< int32_t > >( numerator.uvw[2], level ) )
{
velocityIndices.push_back( dof.value() );
}
}
for ( auto dof : FunctionIterator< vertexdof::VertexDoFFunction< PetscInt > >( numerator.p, level ) )
for ( auto dof : FunctionIterator< vertexdof::VertexDoFFunction< int32_t > >( numerator.p, level ) )
{
pressureIndices.push_back( dof.value() );
}
}
void gatherIndices( std::vector< PetscInt >& velocityIndices,
std::vector< PetscInt >& pressureIndices,
const PrimitiveStorage& storage,
uint_t level,
const P2P1TaylorHoodFunction< PetscInt >& numerator )
void gatherIndices( std::vector< int32_t >& velocityIndices,
std::vector< int32_t >& pressureIndices,
const PrimitiveStorage& storage,
uint_t level,
const P2P1TaylorHoodFunction< int32_t >& numerator )
{
for ( auto dof :
FunctionIterator< vertexdof::VertexDoFFunction< PetscInt > >( numerator.uvw[0].getVertexDoFFunction(), level ) )
FunctionIterator< vertexdof::VertexDoFFunction< int32_t > >( numerator.uvw[0].getVertexDoFFunction(), level ) )
{
velocityIndices.push_back( dof.value() );
}
for ( auto dof : FunctionIterator< EdgeDoFFunction< PetscInt > >( numerator.uvw[0].getEdgeDoFFunction(), level ) )
for ( auto dof : FunctionIterator< EdgeDoFFunction< int32_t > >( numerator.uvw[0].getEdgeDoFFunction(), level ) )
{
velocityIndices.push_back( dof.value() );
}
for ( auto dof :
FunctionIterator< vertexdof::VertexDoFFunction< PetscInt > >( numerator.uvw[1].getVertexDoFFunction(), level ) )
FunctionIterator< vertexdof::VertexDoFFunction< int32_t > >( numerator.uvw[1].getVertexDoFFunction(), level ) )
{
velocityIndices.push_back( dof.value() );
}
for ( auto dof : FunctionIterator< EdgeDoFFunction< PetscInt > >( numerator.uvw[1].getEdgeDoFFunction(), level ) )
for ( auto dof : FunctionIterator< EdgeDoFFunction< int32_t > >( numerator.uvw[1].getEdgeDoFFunction(), level ) )
{
velocityIndices.push_back( dof.value() );
}
if ( storage.hasGlobalCells() )
{
for ( auto dof :
FunctionIterator< vertexdof::VertexDoFFunction< PetscInt > >( numerator.uvw[2].getVertexDoFFunction(), level ) )
FunctionIterator< vertexdof::VertexDoFFunction< int32_t > >( numerator.uvw[2].getVertexDoFFunction(), level ) )
{
velocityIndices.push_back( dof.value() );
}
for ( auto dof : FunctionIterator< EdgeDoFFunction< PetscInt > >( numerator.uvw[2].getEdgeDoFFunction(), level ) )
for ( auto dof : FunctionIterator< EdgeDoFFunction< int32_t > >( numerator.uvw[2].getEdgeDoFFunction(), level ) )
{
velocityIndices.push_back( dof.value() );
}
}
for ( auto dof : FunctionIterator< vertexdof::VertexDoFFunction< PetscInt > >( numerator.p, level ) )
for ( auto dof : FunctionIterator< vertexdof::VertexDoFFunction< int32_t > >( numerator.p, level ) )
{
pressureIndices.push_back( dof.value() );
}
}
bool p2p1StokesPetscApplyTest( const uint_t& level, const std::string& meshFile, const DoFType& location, const real_t& eps )
bool p2p1StokesGinkgoApplyTest( std::shared_ptr< const gko::Executor > exec,
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 );
......@@ -128,14 +130,14 @@ bool p2p1StokesPetscApplyTest( const uint_t& level, const std::string& meshFile,
loadbalancing::roundRobin( setupStorage );
std::shared_ptr< hyteg::PrimitiveStorage > storage = std::make_shared< hyteg::PrimitiveStorage >( setupStorage );
writeDomainPartitioningVTK( storage, "../../output", "P2P1StokesPetscApplyTestDomain" );
writeDomainPartitioningVTK( storage, "../../output", "P2P1StokesGinkgoApplyTestDomain" );
P2P1TaylorHoodFunction< real_t > src( "src", storage, level, level );
P2P1TaylorHoodFunction< real_t > hhgDst( "hhgDst", storage, level, level );
P2P1TaylorHoodFunction< real_t > petscDst( "petscDst", storage, level, level );
P2P1TaylorHoodFunction< real_t > ginkgoDst( "ginkgoDst", storage, level, level );
P2P1TaylorHoodFunction< real_t > err( "error", storage, level, level );
P2P1TaylorHoodFunction< real_t > ones( "ones", storage, level, level );
P2P1TaylorHoodFunction< PetscInt > numerator( "numerator", storage, level, level );
P2P1TaylorHoodFunction< int32_t > 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; };
......@@ -148,7 +150,7 @@ bool p2p1StokesPetscApplyTest( const uint_t& level, const std::string& meshFile,
src.interpolate( srcFunction, level, hyteg::All );
hhgDst.interpolate( rand, level, location );
petscDst.interpolate( rand, level, location );
ginkgoDst.interpolate( rand, level, location );
ones.interpolate( one, level, location );
P2P1TaylorHoodStokesOperator L( storage, level, level );
......@@ -163,11 +165,10 @@ bool p2p1StokesPetscApplyTest( const uint_t& level, const std::string& meshFile,
// HyTeG apply
L.apply( src, hhgDst, level, location );
// PETSc apply
// Ginkgo apply
using mtx = gko::distributed::Matrix< real_t, int32_t >;
using vec = gko::distributed::Vector< real_t, int32_t >;
auto exec = gko::ReferenceExecutor::create();
auto comm = gko::mpi::communicator::create_world();
auto part = gko::share( gko::distributed::Partition< int32_t >::build_uniformly( exec, 1, localDoFs ) );
......@@ -183,7 +184,7 @@ bool p2p1StokesPetscApplyTest( const uint_t& level, const std::string& meshFile,
level,
All );
hyteg::petsc::createVectorFromFunction(
petscDst,
ginkgoDst,
numerator,
std::make_shared< GinkgoVectorProxy >( gko::lend( dstGkoVec ), gko::dim< 2 >{ localDoFs, 1 }, part ),
level,
......@@ -224,14 +225,14 @@ bool p2p1StokesPetscApplyTest( const uint_t& level, const std::string& meshFile,
b_gkoMtx->apply(b_srcGkoVec.get(), b_dstGkoVec.get());
hyteg::petsc::createFunctionFromVector(
petscDst,
ginkgoDst,
numerator,
std::make_shared< GinkgoVectorProxy >( gko::lend( dstGkoVec ), gko::dim< 2 >{ localDoFs, 1 }, part ),
level,
All );
// compare
err.assign( { 1.0, -1.0 }, { hhgDst, petscDst }, level, location );
err.assign( { 1.0, -1.0 }, { hhgDst, ginkgoDst }, level, location );
auto maxError = err.uvw[0].getMaxMagnitude( level );
maxError = std::max( maxError, err.uvw[1].getMaxMagnitude( level ) );
if ( err.uvw.getDimension() == 3 )
......@@ -244,10 +245,10 @@ bool p2p1StokesPetscApplyTest( const uint_t& level, const std::string& meshFile,
if ( writeVTK )
{
VTKOutput vtkOutput( "../../output", "P2P1StokesPetscApplyTest", storage );
VTKOutput vtkOutput( "../../output", "P2P1StokesGinkgoApplyTest", storage );
vtkOutput.add( src.uvw );
vtkOutput.add( hhgDst.uvw );
vtkOutput.add( petscDst.uvw );
vtkOutput.add( ginkgoDst.uvw );
vtkOutput.add( err.uvw );
vtkOutput.write( level, 0 );
}
......@@ -269,17 +270,32 @@ int main( int argc, char* argv[] )
{
walberla::MPIManager::instance()->initializeMPI( &argc, &argv );
walberla::MPIManager::instance()->useWorldComm();
hyteg::PETScManager petscManager( &argc, &argv );
using namespace hyteg;
bool succeeded = true;
succeeded &= hyteg::p2p1StokesPetscApplyTest( 3, "../../data/meshes/quad_4el.msh", hyteg::All, 8.7e-15 );
succeeded &= hyteg::p2p1StokesPetscApplyTest( 3, "../../data/meshes/annulus_coarse.msh", hyteg::All, 2.6e-13 );
succeeded &= hyteg::p2p1StokesPetscApplyTest( 3, "../../data/meshes/3D/tet_1el.msh", hyteg::All, 1.0e-16 );
succeeded &= hyteg::p2p1StokesPetscApplyTest( 2, "../../data/meshes/3D/pyramid_2el.msh", hyteg::All, 7.3e-16 );
succeeded &= hyteg::p2p1StokesPetscApplyTest( 2, "../../data/meshes/3D/pyramid_4el.msh", hyteg::All, 1.4e-15 );
succeeded &= hyteg::p2p1StokesPetscApplyTest( 2, "../../data/meshes/3D/regular_octahedron_8el.msh", hyteg::All, 4.0e-15 );
succeeded &= hyteg::p2p1StokesPetscApplyTest( 2, "../../data/meshes/3D/cube_24el.msh", hyteg::All, 3.5e-15 );
for ( auto tag : tag_list )
{
auto exec = get_executor( tag );
if ( exec )
{
try
{
succeeded &= p2p1StokesGinkgoApplyTest( exec, 3, "../../data/meshes/quad_4el.msh", hyteg::All, 8.7e-15 );
succeeded &= p2p1StokesGinkgoApplyTest( exec, 3, "../../data/meshes/annulus_coarse.msh", hyteg::All, 2.6e-13 );
succeeded &= p2p1StokesGinkgoApplyTest( exec, 3, "../../data/meshes/3D/tet_1el.msh", hyteg::All, 1.0e-16 );
succeeded &= p2p1StokesGinkgoApplyTest( exec, 2, "../../data/meshes/3D/pyramid_2el.msh", hyteg::All, 7.3e-16 );
succeeded &= p2p1StokesGinkgoApplyTest( exec, 2, "../../data/meshes/3D/pyramid_4el.msh", hyteg::All, 1.4e-15 );
succeeded &=
p2p1StokesGinkgoApplyTest( exec, 2, "../../data/meshes/3D/regular_octahedron_8el.msh", hyteg::All, 4.0e-15 );
succeeded &= p2p1StokesGinkgoApplyTest( exec, 2, "../../data/meshes/3D/cube_24el.msh", hyteg::All, 3.5e-15 );
} catch ( const gko::NotImplemented& e )
{
std::cout << e.what() << " for executor " << get_executor_name( exec ) << std::endl;
}
}
}
WALBERLA_CHECK( succeeded, "One of the tests failed" )
......
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