Added small test to check that PETSc and HyTeG apply() do the same thing for the DG operator.

src/hyteg/dgfunctionspace/DGDiffusionForm_Example.hpp

@@ -35,6 +35,14 @@ namespace dg {
 
 using walberla::real_c;
 
+/// \brief Reference implementation of a DG form for the Laplacian.
+///
+/// Refer to
+///     Béatrice Rivière (2008).
+///     Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations: Theory and Implementation
+/// for details.
+///
+/// Note: if uncertain, try setting beta_0 to 1 for 2D and 0.5 for 3D applications.
 class DGDiffusionForm_Example : public DGForm
 {
  public:

src/hyteg/dgfunctionspace/DGFunction.cpp

@@ -758,6 +758,75 @@ void DGFunction< ValueType >::fromVector( const DGFunction< idx_t >&
    WALBERLA_UNUSED( flag );
 }
 
+template < typename ValueType >
+ValueType DGFunction< ValueType >::getMaxMagnitude( uint_t level, bool mpiReduce ) const
+{
+   auto localMax = ValueType( 0.0 );
+
+   if ( storage_->hasGlobalCells() )
+   {
+      for ( auto& it : this->getStorage()->getCells() )
+      {
+         const auto cellID = it.first;
+         const auto cell   = *it.second;
+
+         const auto degree  = polynomialDegree( cellID );
+         const auto numDofs = basis()->numDoFsPerElement( 3, uint_c( degree ) );
+
+         auto       dofs      = volumeDoFFunction()->dofMemory( cellID, level );
+         const auto memLayout = volumeDoFFunction()->memoryLayout();
+
+         for ( auto cellType : celldof::allCellTypes )
+         {
+            for ( const auto& idxIt : celldof::macrocell::Iterator( level, cellType ) )
+            {
+               for ( uint_t i = 0; i < numDofs; i++ )
+               {
+                  const auto val = abs( dofs[volumedofspace::indexing::index(
+                      idxIt.x(), idxIt.y(), idxIt.z(), cellType, i, numDofs, level, memLayout )] );
+                  localMax       = std::max( localMax, val );
+               }
+            }
+         }
+      }
+   }
+   else
+   {
+      for ( auto& it : this->getStorage()->getFaces() )
+      {
+         const auto faceID = it.first;
+         const auto face   = *it.second;
+
+         const auto degree  = polynomialDegree( faceID );
+         const auto numDofs = basis()->numDoFsPerElement( 2, uint_c( degree ) );
+
+         auto       dofs      = volumeDoFFunction()->dofMemory( faceID, level );
+         const auto memLayout = volumeDoFFunction()->memoryLayout();
+
+         for ( auto faceType : facedof::allFaceTypes )
+         {
+            for ( const auto& idxIt : facedof::macroface::Iterator( level, faceType ) )
+            {
+               for ( uint_t i = 0; i < numDofs; i++ )
+               {
+                  const auto val = abs(
+                      dofs[volumedofspace::indexing::index( idxIt.x(), idxIt.y(), faceType, i, numDofs, level, memLayout )] );
+                  localMax = std::max( localMax, val );
+               }
+            }
+         }
+      }
+   }
+
+   ValueType globalMax = localMax;
+   if ( mpiReduce )
+   {
+      globalMax = walberla::mpi::allReduce( localMax, walberla::mpi::MAX );
+   }
+
+   return globalMax;
+}
+
 /// explicit instantiation
 template class DGFunction< double >;
 template class DGFunction< float >;

src/hyteg/dgfunctionspace/DGFunction.hpp

@@ -337,6 +337,12 @@ class DGFunction final : public Function< DGFunction< ValueType > >
       WALBERLA_ABORT( "DGFunction::swap() not implemented." )
    }
 
+   /// \brief Returns the max absolute DoF.
+   ///
+   /// \param level     refinement level
+   /// \param mpiReduce if true, reduces over all processes (global max magnitude), if false returns the process local value
+   ValueType getMaxMagnitude( uint_t level, bool mpiReduce = true ) const;
+
  private:
    using Function< DGFunction< ValueType > >::communicators_;
    using Function< DGFunction< ValueType > >::additiveCommunicators_;

tests/hyteg/CMakeLists.txt

@@ -798,6 +798,9 @@ if (HYTEG_BUILD_WITH_PETSC)
 
     waLBerla_compile_test(FILES dg/DGInterpolateEvaluateTest.cpp DEPENDS hyteg core)
     waLBerla_execute_test(NAME DGInterpolateEvaluateTest)
+
+    waLBerla_compile_test(FILES dg/DGPetscApplyTest.cpp DEPENDS hyteg core)
+    waLBerla_execute_test(NAME DGPetscApplyTest)
 endif()
 
 ## Basic Smoothing and Solving

tests/hyteg/dg/DGPetscApplyTest.cpp

+/*
+* Copyright (c) 2017-2022 Nils Kohl.
+*
+* 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/DataTypes.h"
+#include "core/math/Random.h"
+#include "core/mpi/MPIManager.h"
+
+#include "hyteg/dataexport/VTKOutput.hpp"
+#include "hyteg/dgfunctionspace/DGBasisLinearLagrange_Example.hpp"
+#include "hyteg/dgfunctionspace/DGDiffusionForm_Example.hpp"
+#include "hyteg/dgfunctionspace/DGFunction.hpp"
+#include "hyteg/dgfunctionspace/DGMassForm_Example.hpp"
+#include "hyteg/dgfunctionspace/DGOperator.hpp"
+#include "hyteg/functions/FunctionTraits.hpp"
+#include "hyteg/mesh/MeshInfo.hpp"
+#include "hyteg/petsc/PETScCGSolver.hpp"
+#include "hyteg/petsc/PETScManager.hpp"
+#include "hyteg/petsc/PETScSparseMatrix.hpp"
+#include "hyteg/petsc/PETScVector.hpp"
+#include "hyteg/primitivestorage/PrimitiveStorage.hpp"
+#include "hyteg/primitivestorage/SetupPrimitiveStorage.hpp"
+
+using walberla::real_t;
+
+namespace hyteg {
+
+void dgPetscApplyTest( uint_t level, const MeshInfo& meshInfo, real_t eps )
+{
+   using namespace dg;
+
+   SetupPrimitiveStorage setupStorage( meshInfo, walberla::uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) );
+   setupStorage.setMeshBoundaryFlagsOnBoundary( 1, 0, true );
+   auto storage = std::make_shared< PrimitiveStorage >( setupStorage, 1 );
+
+   auto basis    = std::make_shared< DGBasisLinearLagrange_Example >();
+   auto diffForm = std::make_shared< DGDiffusionForm_Example >( ( storage->hasGlobalCells() ? 0.5 : 1 ) );
+   auto massForm = std::make_shared< DGMassForm_Example >();
+
+   DGFunction< real_t > src( "src", storage, level, level, basis, 1 );
+   DGFunction< real_t > tmp( "tmp", storage, level, level, basis, 1 );
+   DGFunction< real_t > hytegDst( "hytegDst", storage, level, level, basis, 1 );
+   DGFunction< real_t > petscDst( "petscDst", storage, level, level, basis, 1 );
+   DGFunction< real_t > err( "error", storage, level, level, basis, 1 );
+   DGFunction< idx_t >  numerator( "numerator", storage, level, level, basis, 1 );
+
+   DGOperator L( storage, level, level, diffForm );
+   DGOperator M( storage, level, level, massForm );
+
+   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] );
+   };
+
+   // interpolation of some function into the src vector
+   tmp.evaluateLinearFunctional( srcFunction, level );
+   PETScCGSolver< DGOperator > solverM( storage, level, numerator );
+   solverM.solve( M, src, tmp, level );
+
+   numerator.enumerate( level );
+
+   const uint_t globalDoFs = src.getNumberOfGlobalDoFs( level );
+   WALBERLA_LOG_INFO_ON_ROOT( "Global DoFs: " << globalDoFs );
+
+   // HyTeG apply
+   L.apply( src, hytegDst, level, All, Replace );
+
+   // PETSc apply
+   PETScVector< real_t, DGFunction > srcPetscVec;
+   PETScVector< real_t, DGFunction > dstPetscVec;
+   PETScSparseMatrix< DGOperator >   petscMatrix;
+
+   srcPetscVec.createVectorFromFunction( src, numerator, level );
+   dstPetscVec.createVectorFromFunction( petscDst, numerator, level );
+   petscMatrix.createMatrixFromOperator( L, level, numerator );
+
+   WALBERLA_CHECK( petscMatrix.isSymmetric() );
+
+   MatMult( petscMatrix.get(), srcPetscVec.get(), dstPetscVec.get() );
+
+   dstPetscVec.createFunctionFromVector( petscDst, numerator, level );
+
+   // compare
+   err.assign( { 1.0, -1.0 }, { hytegDst, petscDst }, level );
+   auto maxMag = err.getMaxMagnitude( level );
+
+   WALBERLA_LOG_INFO_ON_ROOT( "Error max mag = " << maxMag );
+
+   // VTK
+   //   VTKOutput vtkOutput( "../../output", "DGPetscApplyTest", storage );
+   //   vtkOutput.add( src );
+   //   vtkOutput.add( hytegDst );
+   //   vtkOutput.add( petscDst );
+   //   vtkOutput.add( err );
+   //   vtkOutput.write( level, 0 );
+
+   WALBERLA_CHECK_LESS( maxMag, eps );
+}
+
+} // namespace hyteg
+
+int main( int argc, char* argv[] )
+{
+   walberla::MPIManager::instance()->initializeMPI( &argc, &argv );
+   walberla::MPIManager::instance()->useWorldComm();
+   hyteg::PETScManager petscManager( &argc, &argv );
+
+   hyteg::dgPetscApplyTest( 3, hyteg::MeshInfo::fromGmshFile( "../../data/meshes/quad_4el.msh" ), 4.0e-15 );
+   hyteg::dgPetscApplyTest( 3, hyteg::MeshInfo::fromGmshFile( "../../data/meshes/annulus_coarse.msh" ), 6.0e-14 );
+   hyteg::dgPetscApplyTest( 3, hyteg::MeshInfo::fromGmshFile( "../../data/meshes/3D/tet_1el.msh" ), 2.0e-17 );
+   hyteg::dgPetscApplyTest( 3, hyteg::MeshInfo::fromGmshFile( "../../data/meshes/3D/pyramid_2el.msh" ), 5.0e-16 );
+   hyteg::dgPetscApplyTest( 3, hyteg::MeshInfo::fromGmshFile( "../../data/meshes/3D/pyramid_4el.msh" ), 1.0e-15 );
+   hyteg::dgPetscApplyTest( 3, hyteg::MeshInfo::fromGmshFile( "../../data/meshes/3D/regular_octahedron_8el.msh" ), 1.0e-15 );
+
+   return EXIT_SUCCESS;
+}