[wip] add ginkgo blocked solver

src/hyteg/ginkgo/GinkgoBlockSolver.hpp

+#pragma once
+
+#include <ginkgo/core/log/convergence.hpp>
+#include <ginkgo/core/log/convergence_stream.hpp>
+#include <ginkgo/core/matrix/block_matrix.hpp>
+#include <ginkgo/core/solver/bicgstab.hpp>
+#include <ginkgo/core/solver/gmres.hpp>
+#include <ginkgo/core/stop/iteration.hpp>
+#include <ginkgo/core/stop/residual_norm.hpp>
+#include <hyteg/composites/petsc/P2P1TaylorHoodPetsc.hpp>
+
+#include "hyteg/ginkgo/GinkgoCommunication.hpp"
+#include "hyteg/ginkgo/GinkgoDirichletHandling.hpp"
+#include "hyteg/ginkgo/GinkgoSparseMatrixProxy.hpp"
+#include "hyteg/ginkgo/GinkgoVectorProxy.hpp"
+#include "hyteg/solvers/Solver.hpp"
+
+namespace hyteg {
+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< int32_t > >( numerator.uvw[0], level ) )
+   {
+      velocityIndices.push_back( dof.value() );
+   }
+
+   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< int32_t > >( numerator.uvw[2], level ) )
+      {
+         velocityIndices.push_back( dof.value() );
+      }
+   }
+   for ( auto dof : FunctionIterator< vertexdof::VertexDoFFunction< int32_t > >( numerator.p, level ) )
+   {
+      pressureIndices.push_back( dof.value() );
+   }
+}
+
+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< int32_t > >( numerator.uvw[0].getVertexDoFFunction(), level ) )
+   {
+      velocityIndices.push_back( dof.value() );
+   }
+   for ( auto dof : FunctionIterator< EdgeDoFFunction< int32_t > >( numerator.uvw[0].getEdgeDoFFunction(), level ) )
+   {
+      velocityIndices.push_back( dof.value() );
+   }
+   for ( auto dof :
+         FunctionIterator< vertexdof::VertexDoFFunction< int32_t > >( numerator.uvw[1].getVertexDoFFunction(), level ) )
+   {
+      velocityIndices.push_back( dof.value() );
+   }
+   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< int32_t > >( numerator.uvw[2].getVertexDoFFunction(), level ) )
+      {
+         velocityIndices.push_back( dof.value() );
+      }
+      for ( auto dof : FunctionIterator< EdgeDoFFunction< int32_t > >( numerator.uvw[2].getEdgeDoFFunction(), level ) )
+      {
+         velocityIndices.push_back( dof.value() );
+      }
+   }
+   for ( auto dof : FunctionIterator< vertexdof::VertexDoFFunction< int32_t > >( numerator.p, level ) )
+   {
+      pressureIndices.push_back( dof.value() );
+   }
+}
+
+
+template < typename OperatorType >
+class GinkgoBlockSolver : public Solver< OperatorType >
+{
+   using BlockPreconditioner_T = typename OperatorType::BlockPreconditioner_T;
+
+ public:
+   using FunctionType = typename OperatorType::srcType;
+   using valueType    = typename FunctionType::valueType;
+
+   using mtx   = gko::distributed::Matrix< valueType, int32_t >;
+   using csr   = gko::matrix::Csr< valueType, int32_t >;
+   using vec   = gko::distributed::Vector< valueType, int32_t >;
+   using dense = gko::matrix::Dense< valueType >;
+
+   explicit GinkgoBlockSolver( const std::shared_ptr< PrimitiveStorage >& storage,
+                               const uint_t&                              level,
+                               std::shared_ptr< const gko::Executor >     exec )
+   : storage_( storage )
+   , exec_( exec )
+   , comm_( gko::mpi::communicator::create( storage->getSplitCommunicatorByPrimitiveDistribution() ) )
+   , num_( "numerator", storage, level, level )
+   {}
+
+   void solve( const OperatorType&                   A,
+               const typename OperatorType::srcType& x,
+               const typename OperatorType::dstType& b,
+               const walberla::uint_t                level ) override
+   {
+      const auto num_local_dofs  = numberOfLocalDoFs< typename FunctionType::Tag >( *storage_, level );
+      const auto num_global_dofs = numberOfGlobalDoFs< typename FunctionType::Tag >( *storage_, level, comm_->get() );
+
+      auto [start, end] = local_range( num_local_dofs, comm_ );
+      auto part         = gko::share( gko::distributed::Partition<>::build_from_local_range( exec_, start, end, comm_ ) );
+
+      // maybe called in parallel, thus need to keep it for empty processes
+      num_.copyBoundaryConditionFromFunction( x );
+      num_.enumerate( level );
+
+      std::vector< int32_t > bcIndices;
+      hyteg::petsc::applyDirichletBC( num_, bcIndices, level );
+      std::sort( std::begin( bcIndices ), std::end( bcIndices ) );
+
+      auto x_gko = vec::create( exec_, comm_, part, gko::dim< 2 >{ num_global_dofs, 1 }, gko::dim< 2 >{ num_local_dofs, 1 } );
+      auto b_gko = vec::create( exec_, comm_, part, gko::dim< 2 >{ num_global_dofs, 1 }, gko::dim< 2 >{ num_local_dofs, 1 } );
+
+      hyteg::petsc::createVectorFromFunction(
+          x, num_, std::make_shared< GinkgoVectorProxy >( x_gko.get(), gko::dim< 2 >{ num_global_dofs, 1 }, part ), level, All );
+
+      hyteg::petsc::createVectorFromFunction(
+          b, num_, std::make_shared< GinkgoVectorProxy >( b_gko.get(), gko::dim< 2 >{ num_global_dofs, 1 }, part ), level, All );
+
+      monolithic_matrix_ = gko::share( mtx::create( exec_, comm_ ) );
+
+      auto proxy = std::make_shared< GinkgoSparseMatrixProxy< mtx > >(
+          monolithic_matrix_.get(), gko::dim< 2 >{ num_global_dofs, num_global_dofs }, part );
+      hyteg::petsc::createMatrix( A, num_, num_, proxy, level, All );
+      proxy->finalize();
+
+      auto dir_handler = std::make_unique< ZeroRowsDirichletHandler >( bcIndices, b_gko.get(), monolithic_matrix_, true );
+      dir_handler->update_matrix();
+
+      std::vector< int32_t > vIndices;
+      std::vector< int32_t > pIndices;
+
+      gatherIndices( vIndices, pIndices, *storage_, level, num_ );
+
+      auto [p_v_min, p_v_max] = std::minmax_element( std::begin( vIndices ), std::end( vIndices ) );
+      auto [p_p_min, p_p_max] = std::minmax_element( std::begin( pIndices ), std::end( pIndices ) );
+      gko::span v_span( *p_v_min, *p_v_max + 1 );
+      gko::span p_span( *p_p_min, *p_p_max + 1 );
+      if ( p_span.begin != v_span.end )
+      {
+         WALBERLA_ABORT("Indices are NOT blocked: v" << v_span << ", p" << p_span)
+      }
+
+      auto b_mtx_v  = gko::share( monolithic_matrix_->get_local_diag()->create_submatrix( v_span, v_span ) );
+      auto b_mtx_vp = gko::share( monolithic_matrix_->get_local_diag()->create_submatrix( v_span, p_span ) );
+      auto b_mtx_pv = gko::share( monolithic_matrix_->get_local_diag()->create_submatrix( p_span, v_span ) );
+      auto b_mtx_p  = gko::share( monolithic_matrix_->get_local_diag()->create_submatrix( p_span, p_span ) );
+      matrix_ = gko::share(gko::matrix::BlockMatrix::create(exec_, {num_global_dofs, num_global_dofs}, {{b_mtx_v, b_mtx_vp}, {b_mtx_pv, b_mtx_p}}));
+
+//      auto b_vec_b0 = gko::share( b_gko->get_local()->create_submatrix( v_span, 0 ) );
+//      auto b_vec_b1 = gko::share( b_gko->get_local()->create_submatrix( p_span, 0 ) );
+//      auto b_vec_x0 = gko::share( x_gko->get_local()->create_submatrix( v_span, 0 ) );
+//      auto b_vec_x1 = gko::share( x_gko->get_local()->create_submatrix( p_span, 0 ) );
+//
+//      auto b_b_gko = gko::matrix::BlockMatrix::create(exec_, {num_global_dofs, num_global_dofs}, {{b_vec_b0}, {b_vec_b1}});
+//      auto b_x_gko = gko::matrix::BlockMatrix::create(exec_, {num_global_dofs, num_global_dofs}, {{b_vec_x0}, {b_vec_x1}});
+
+      solver_ = gko::solver::Bicgstab< valueType >::build()
+                    .with_criteria( gko::share( gko::stop::ResidualNorm<>::build()
+                                                    .with_baseline( gko::stop::mode::initial_resnorm )
+                                                    .with_reduction_factor( 1e-30 )
+                                                    .on( exec_ ) ),
+                                    gko::share( gko::stop::ResidualNorm<>::build()
+                                                    .with_baseline( gko::stop::mode::absolute )
+                                                    .with_reduction_factor( 1e-12 )
+                                                    .on( exec_ ) ),
+                                    gko::share( gko::stop::Iteration::build().with_max_iters( 15000 ).on( exec_ ) ) )
+                    .on( exec_ )
+                    ->generate( matrix_ );
+
+      auto log     = gko::share( gko::log::Convergence< valueType >::create( exec_ ) );
+      auto factory = gko::clone( solver_->get_stop_criterion_factory() );
+      factory->add_logger( log );
+      solver_->set_stop_criterion_factory( gko::share( factory ) );
+
+      auto rhs = dir_handler->get_rhs( b_gko.get(), x_gko.get() );
+      auto x0  = dir_handler->get_initial_guess( x_gko.get() );
+      solver_->apply( rhs, x0 );
+      dir_handler->update_solution( x0 );
+
+      hyteg::petsc::createFunctionFromVector(
+          x, num_, std::make_shared< GinkgoVectorProxy >( x0, gko::dim< 2 >{ num_global_dofs, 1 }, part ), level, All );
+
+      {
+         WALBERLA_LOG_INFO_ON_ROOT(
+             "[Ginkgo CG]" << ( !log->has_converged() ? " NOT " : " " ) << "converged after " << log->get_num_iterations()
+                           << " iterations,  residual norm: "
+                           << exec_->copy_val_to_host(
+                                  gko::as< gko::matrix::Dense< valueType > >( log->get_residual_norm() )->get_const_values() ) );
+      }
+   }
+
+ private:
+   std::shared_ptr< PrimitiveStorage >                   storage_;
+   std::shared_ptr< const gko::Executor >                exec_;
+   std::shared_ptr< gko::mpi::communicator >             comm_;
+   std::unique_ptr< gko::solver::Bicgstab< valueType > > solver_;
+
+   typename OperatorType::srcType::template FunctionType< int > num_;
+
+   std::shared_ptr< gko::matrix::BlockMatrix > matrix_;
+   std::shared_ptr< mtx >                      monolithic_matrix_;
+};
+
+} // namespace hyteg

tests/hyteg/CMakeLists.txt

@@ -352,11 +352,18 @@ if( HYTEG_BUILD_WITH_PETSC )
 endif()
 
 if( HYTEG_BUILD_WITH_GINKGO )
+  add_library(precompiled composites/p2p1stokes.cpp)
+  target_link_modules(precompiled hyteg core)
+  target_link_libraries(precompiled ${SERVICE_LIBS})
+
   waLBerla_compile_test(FILES P1/P1GinkgoSolveTest.cpp DEPENDS hyteg core)
   waLBerla_execute_test(NAME P1GinkgoSolveTest1 COMMAND $<TARGET_FILE:P1GinkgoSolveTest> PROCESSES 1)
 
   waLBerla_compile_test(FILES composites/P2P1StokesGinkgoApplyTest.cpp DEPENDS hyteg core)
   waLBerla_execute_test(NAME P2P1StokesGinkgoApplyTest1 COMMAND $<TARGET_FILE:P2P1StokesGinkgoApplyTest> PROCESSES 1)
+
+  waLBerla_compile_test(FILES composites/P2P1Stokes3DGinkgoSolveTest.cpp DEPENDS hyteg core)
+  waLBerla_execute_test(NAME P2P1Stokes3DGinkgoSolveTest1 COMMAND $<TARGET_FILE:P2P1Stokes3DGinkgoSolveTest> PROCESSES 1)
 endif()
 
 

tests/hyteg/composites/P2P1Stokes3DGinkgoSolveTest.cpp

+/*
+ * Copyright (c) 2017-2019 Dominik Thoennes.
+ *
+ * This file is part of HyTeG
+ * (see https://i10git.cs.fau.de/hyteg/hyteg).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * 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 "core/Environment.h"
+#include "core/logging/Logging.h"
+#include "core/math/Random.h"
+#include "core/timing/Timer.h"
+
+#include "hyteg/dataexport/VTKOutput.hpp"
+#include "hyteg/mesh/MeshInfo.hpp"
+#include "hyteg/misc/ExactStencilWeights.hpp"
+#include "hyteg/p1functionspace/P1ConstantOperator.hpp"
+#include "hyteg/p1functionspace/P1Function.hpp"
+#include "hyteg/ginkgo/GinkgoBlockSolver.hpp"
+#include "hyteg/petsc/PETScLUSolver.hpp"
+#include "hyteg/petsc/PETScMinResSolver.hpp"
+#include "hyteg/petsc/PETScBlockPreconditionedStokesSolver.hpp"
+#include "hyteg/petsc/PETScManager.hpp"
+#include "hyteg/petsc/PETScVersion.hpp"
+#include "hyteg/petsc/PETScExportLinearSystem.hpp"
+#include "hyteg/ginkgo/GinkgoSparseMatrixProxy.hpp"
+#include "hyteg/ginkgo/GinkgoUtilities.hpp"
+#include "hyteg/ginkgo/GinkgoVectorProxy.hpp"
+#include "hyteg/primitivestorage/SetupPrimitiveStorage.hpp"
+#include "hyteg/primitivestorage/Visualization.hpp"
+#include "hyteg/primitivestorage/loadbalancing/SimpleBalancer.hpp"
+#include "hyteg/composites/P2P1TaylorHoodFunction.hpp"
+#include "hyteg/composites/P2P1TaylorHoodStokesOperator.hpp"
+
+#ifndef HYTEG_BUILD_WITH_GINKGO
+  WALBERLA_ABORT( "This test only works with Ginkgo enabled. Please enable it via -DHYTEG_BUILD_WITH_GINKGO=ON" )
+#endif
+
+using walberla::real_t;
+using walberla::uint_c;
+using walberla::uint_t;
+
+namespace hyteg {
+
+void petscSolveTest(const uint_t & level, const MeshInfo & meshInfo,  const real_t & resEps, const real_t & errEpsUSum, const real_t & errEpsP )
+{
+  SetupPrimitiveStorage setupStorage( meshInfo, uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) );
+
+  setupStorage.setMeshBoundaryFlagsOnBoundary( 1, 0, true );
+
+  hyteg::loadbalancing::roundRobin( setupStorage );
+
+  std::shared_ptr< PrimitiveStorage > storage = std::make_shared< PrimitiveStorage >( setupStorage );
+  writeDomainPartitioningVTK( storage, "../../output", "P2P1Stokes3DPetscSolve_Domain" );
+
+  hyteg::P2P1TaylorHoodFunction< real_t >                      x( "x", storage, level, level );
+  hyteg::P2P1TaylorHoodFunction< real_t >                      x_exact( "x_exact", storage, level, level );
+  hyteg::P2P1TaylorHoodFunction< real_t >                      b( "b", storage, level, level );
+  hyteg::P2P1TaylorHoodFunction< real_t >                      btmp( "btmp", storage, level, level );
+  hyteg::P2P1TaylorHoodFunction< real_t >                      err( "err", storage, level, level );
+  hyteg::P2P1TaylorHoodFunction< real_t >                      residuum( "res", storage, level, level );
+  hyteg::P2P1TaylorHoodFunction< real_t >                      nullspace( "nullspace", storage, level, level );
+  hyteg::P2P1TaylorHoodFunction< PetscInt >                    numerator( "numerator", storage, level, level );
+
+  numerator.enumerate( level );
+
+  hyteg::P2P1TaylorHoodStokesOperator A( storage, level, level );
+
+  std::function< real_t( const hyteg::Point3D& ) > exactU = []( const hyteg::Point3D& xx ) { return real_c(20) * xx[0] * xx[1] * xx[1] * xx[1]; };
+  std::function< real_t( const hyteg::Point3D& ) > exactV = []( const hyteg::Point3D& xx ) { return real_c(5) * xx[0] * xx[0] * xx[0] * xx[0] - real_c(5) * xx[1] * xx[1] * xx[1] * xx[1]; };
+  std::function< real_t( const hyteg::Point3D& ) > exactW = []( const hyteg::Point3D&    ) { return real_c(0); };
+  std::function< real_t( const hyteg::Point3D& ) > exactP = []( const hyteg::Point3D& xx ) { return real_c(60) * std::pow( xx[0], 2.0 ) * xx[1] - real_c(20) * std::pow( xx[1], 3.0 ); };
+  std::function< real_t( const hyteg::Point3D& ) > zero =   []( const hyteg::Point3D&    ) { return real_c(0); };
+  std::function< real_t( const hyteg::Point3D& ) > ones =   []( const hyteg::Point3D&    ) { return real_c(1); };
+
+  b.uvw.interpolate( { exactU, exactV, exactW }, level, DirichletBoundary );
+  b.p.interpolate( zero, level, All );
+
+  x.uvw.interpolate( { exactU, exactV, exactW }, level, DirichletBoundary );
+
+  x_exact.uvw.interpolate( { exactU, exactV, exactW }, level );
+  x_exact.p.interpolate( exactP, level );
+
+  hyteg::vertexdof::projectMean( x_exact.p, level );
+
+  nullspace.p.interpolate( ones, level, All );
+
+  uint_t localDoFs1 = hyteg::numberOfLocalDoFs< P2P1TaylorHoodFunctionTag >( *storage, level );
+  uint_t globalDoFs1 = hyteg::numberOfGlobalDoFs< P2P1TaylorHoodFunctionTag >( *storage, level );
+  uint_t globalDoFsvelocity = 3 * hyteg::numberOfGlobalDoFs< P2FunctionTag >( *storage, level );
+
+  WALBERLA_LOG_INFO( "localDoFs: " << localDoFs1 << " globalDoFs: " << globalDoFs1
+                                   << ", global velocity dofs: " << globalDoFsvelocity );
+
+  GinkgoBlockSolver<P2P1TaylorHoodStokesOperator> solver(storage, level, gko::ReferenceExecutor::create());
+
+  walberla::WcTimer timer;
+  solver.solve(A, x, b, level);
+  timer.end();
+
+  hyteg::vertexdof::projectMean( x.p, level );
+
+  WALBERLA_LOG_INFO_ON_ROOT( "time was: " << timer.last() );
+  A.apply( x, residuum, level, hyteg::Inner );
+
+  err.assign( {1.0, -1.0}, {x, x_exact}, level );
+
+  real_t discr_l2_err = std::sqrt( err.dotGlobal( err, level ) / (real_t) globalDoFs1 );
+  real_t discr_l2_err_1_u = std::sqrt( err.uvw[0].dotGlobal( err.uvw[0], level ) / (real_t) globalDoFs1 );
+  real_t discr_l2_err_1_v = std::sqrt( err.uvw[1].dotGlobal( err.uvw[1], level ) / (real_t) globalDoFs1 );
+  real_t discr_l2_err_1_w = std::sqrt( err.uvw[2].dotGlobal( err.uvw[2], level ) / (real_t) globalDoFs1 );
+  real_t discr_l2_err_1_p = std::sqrt( err.p.dotGlobal( err.p, level ) / (real_t) globalDoFs1 );
+  real_t residuum_l2_1  = std::sqrt( residuum.dotGlobal( residuum, level ) / (real_t) globalDoFs1 );
+
+  WALBERLA_LOG_INFO_ON_ROOT( "discrete L2 error = " << discr_l2_err );
+  WALBERLA_LOG_INFO_ON_ROOT( "discrete L2 error u = " << discr_l2_err_1_u );
+  WALBERLA_LOG_INFO_ON_ROOT( "discrete L2 error v = " << discr_l2_err_1_v );
+  WALBERLA_LOG_INFO_ON_ROOT( "discrete L2 error w = " << discr_l2_err_1_w );
+  WALBERLA_LOG_INFO_ON_ROOT( "discrete L2 error p = " << discr_l2_err_1_p );
+  WALBERLA_LOG_INFO_ON_ROOT( "residuum 1 = " << residuum_l2_1 );
+
+  WALBERLA_CHECK_LESS( residuum_l2_1, resEps );
+  WALBERLA_CHECK_LESS( discr_l2_err_1_u + discr_l2_err_1_v + discr_l2_err_1_w, errEpsUSum );
+  WALBERLA_CHECK_LESS( discr_l2_err_1_p, errEpsP);
+
+  auto tt = storage->getTimingTree()->getReduced().getCopyWithRemainder();
+}
+
+}
+
+using namespace hyteg;
+
+int main( int argc, char* argv[] )
+{
+  walberla::Environment walberlaEnv( argc, argv );
+  walberla::MPIManager::instance()->useWorldComm();
+
+  PETScManager petscManager( &argc, &argv );
+
+  printPETScVersionNumberString();
+
+  auto level = argc > 2 ? std::stoi(argv[2]) : 0;
+
+  petscSolveTest( level, hyteg::MeshInfo::fromGmshFile( "../../data/meshes/3D/cube_center_at_origin_24el.msh" ), 2.9e-12, 0.021, 0.33 );
+
+  return EXIT_SUCCESS;
+}

tests/hyteg/composites/P2P1StokesGinkgoApplyTest.cpp

@@ -114,6 +114,35 @@ void gatherIndices( std::vector< int32_t >&                  velocityIndices,
    }
 }
 
+bool check_error( std::shared_ptr< PrimitiveStorage >     storage,
+                  const P2P1TaylorHoodFunction< real_t >& a,
+                  const P2P1TaylorHoodFunction< real_t >& b,
+                  const uint                              level,
+                  const real_t&                           eps )
+{ // compare
+   P2P1TaylorHoodFunction< real_t > err( "error", storage, level, level );
+   err.assign( { 1.0, -1.0 }, { a, b }, level, hyteg::All );
+   auto maxError = err.uvw[0].getMaxMagnitude( level );
+   maxError      = std::max( maxError, err.uvw[1].getMaxMagnitude( level ) );
+   if ( err.uvw.getDimension() == 3 )
+   {
+      maxError = std::max( maxError, err.uvw[2].getMaxMagnitude( level ) );
+   }
+   maxError = std::max( maxError, err.p.getMaxMagnitude( level ) );
+
+   WALBERLA_LOG_INFO_ON_ROOT( "Error max Magnitude = " << maxError << " eps: " << eps );
+
+   if ( maxError > eps )
+   {
+      WALBERLA_LOG_INFO_ON_ROOT( "TEST FAILED!" );
+      return false;
+   }
+   else
+   {
+      return true;
+   }
+}
+
 bool p2p1StokesGinkgoApplyTest( std::shared_ptr< const gko::Executor > exec,
                                 const uint_t&                          level,
                                 const std::string&                     meshFile,
@@ -194,6 +223,18 @@ bool p2p1StokesGinkgoApplyTest( std::shared_ptr< const gko::Executor > exec,
    hyteg::petsc::createMatrix( L, numerator, numerator, proxy, level, All );
    proxy->finalize();
 
+   gkoMtx->apply(srcGkoVec.get(), dstGkoVec.get());
+
+   hyteg::petsc::createFunctionFromVector(
+       ginkgoDst,
+       numerator,
+       std::make_shared< GinkgoVectorProxy >( gko::lend( dstGkoVec ), gko::dim< 2 >{ localDoFs, 1 }, part ),
+       level,
+       All );
+
+   auto success = check_error(storage, hhgDst, ginkgoDst, level , eps);
+
+
    std::vector< int32_t > vIndices;
    std::vector< int32_t > pIndices;
 
@@ -222,6 +263,13 @@ bool p2p1StokesGinkgoApplyTest( std::shared_ptr< const gko::Executor > exec,
    auto b_srcGkoVec = gko::matrix::BlockMatrix::create(exec, {localDoFs, localDoFs}, {{b_vec_s0}, {b_vec_s1}});
    auto b_dstGkoVec = gko::matrix::BlockMatrix::create(exec, {localDoFs, localDoFs}, {{b_vec_d0}, {b_vec_d1}});
 
+   hyteg::petsc::createVectorFromFunction(
+       ginkgoDst,
+       numerator,
+       std::make_shared< GinkgoVectorProxy >( gko::lend( dstGkoVec ), gko::dim< 2 >{ localDoFs, 1 }, part ),
+       level,
+       All );
+
    b_gkoMtx->apply(b_srcGkoVec.get(), b_dstGkoVec.get());
 
    hyteg::petsc::createFunctionFromVector(
@@ -232,36 +280,9 @@ bool p2p1StokesGinkgoApplyTest( std::shared_ptr< const gko::Executor > exec,
        All );
 
    // compare
-   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 )
-   {
-      maxError = std::max( maxError, err.uvw[2].getMaxMagnitude( level ) );
-   }
-   maxError = std::max( maxError, err.p.getMaxMagnitude( level ) );
+   success &= check_error(storage, hhgDst, ginkgoDst, level, eps);
 
-   WALBERLA_LOG_INFO_ON_ROOT( "Error max Magnitude = " << maxError << " eps: " << eps );
-
-   if ( writeVTK )
-   {
-      VTKOutput vtkOutput( "../../output", "P2P1StokesGinkgoApplyTest", storage );
-      vtkOutput.add( src.uvw );
-      vtkOutput.add( hhgDst.uvw );
-      vtkOutput.add( ginkgoDst.uvw );
-      vtkOutput.add( err.uvw );
-      vtkOutput.write( level, 0 );
-   }
-
-   if ( maxError > eps )
-   {
-      WALBERLA_LOG_INFO_ON_ROOT( "TEST FAILED!" );
-      return false;
-   }
-   else
-   {
-      return true;
-   }
+   return success;
 }
 
 } // namespace hyteg