Merge branch 'mohr/BlockFunctionToVector-mc' into 'master'

Implements PETSc vector <-> function conversion for Block/GenericFunction

See merge request !447

src/hyteg/dgfunctionspace/DGPetsc.hpp

+/*
+ * Copyright (c) 2021 Marcus Mohr.
+ *
+ * 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/>.
+ */
+#pragma once
+
+#include "hyteg/dgfunctionspace/DGFunction.hpp"
+#include "hyteg/petsc/PETScWrapper.hpp"
+#include "hyteg/sparseassembly/SparseMatrixProxy.hpp"
+#include "hyteg/sparseassembly/VectorProxy.hpp"
+
+#ifdef HYTEG_BUILD_WITH_PETSC
+
+namespace hyteg {
+namespace petsc {
+
+inline void createVectorFromFunction( const DGFunction< PetscReal >&        function,
+                                      const DGFunction< PetscInt >&         numerator,
+                                      const std::shared_ptr< VectorProxy >& vec,
+                                      uint_t                                level,
+                                      DoFType                               flag )
+{
+  WALBERLA_ABORT( "Congrats :( You have detected another unimplemented feature of DGFunction" );
+}
+
+inline void createFunctionFromVector( const DGFunction< PetscReal >&        function,
+                                      const DGFunction< PetscInt >&         numerator,
+                                      const std::shared_ptr< VectorProxy >& vec,
+                                      uint_t                                level,
+                                      DoFType                               flag )
+{
+  WALBERLA_ABORT( "Congrats :( You have detected another unimplemented feature of DGFunction" );
+}
+
+} // namespace petsc
+} // namespace hyteg
+
+#endif

src/hyteg/facedofspace/FaceDoFPetsc.hpp

+/*
+ * Copyright (c) 2021 Marcus Mohr.
+ *
+ * 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/>.
+ */
+#pragma once
+
+#include "hyteg/petsc/PETScWrapper.hpp"
+#include "hyteg/sparseassembly/SparseMatrixProxy.hpp"
+#include "hyteg/sparseassembly/VectorProxy.hpp"
+
+#ifdef HYTEG_BUILD_WITH_PETSC
+
+namespace hyteg {
+namespace petsc {
+
+inline void createVectorFromFunction( const FaceDoFFunction< PetscReal >&        function,
+                                      const FaceDoFFunction< PetscInt >&         numerator,
+                                      const std::shared_ptr< VectorProxy >& vec,
+                                      uint_t                                level,
+                                      DoFType                               flag )
+{
+  WALBERLA_ABORT( "Congrats :( You have detected another unimplemented feature of FaceDoFFunction" );
+}
+
+inline void createFunctionFromVector( const FaceDoFFunction< PetscReal >&        function,
+                                      const FaceDoFFunction< PetscInt >&         numerator,
+                                      const std::shared_ptr< VectorProxy >& vec,
+                                      uint_t                                level,
+                                      DoFType                               flag )
+{
+  WALBERLA_ABORT( "Congrats :( You have detected another unimplemented feature of FaceDoFFunction" );
+}
+
+} // namespace petsc
+} // namespace hyteg
+
+#endif

src/hyteg/functions/BlockFunctionPetsc.hpp

+/*
+ * Copyright (c) 2021 Marcus Mohr.
+ *
+ * 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/>.
+ */
+
+#pragma once
+
+#include "core/DataTypes.h"
+
+#include "hyteg/functions/BlockFunction.hpp"
+#include "hyteg/functions/GenericFunctionPetsc.hpp"
+#include "hyteg/petsc/PETScWrapper.hpp"
+#include "hyteg/sparseassembly/VectorProxy.hpp"
+#include "hyteg/types/flags.hpp"
+
+#ifdef HYTEG_BUILD_WITH_PETSC
+namespace hyteg {
+namespace petsc {
+
+inline void createVectorFromFunction( const BlockFunction< PetscReal >&   function,
+                                      const BlockFunction< PetscInt >&    numerator,
+                                      const std::shared_ptr< VectorProxy >& vec,
+                                      uint_t                                level,
+                                      DoFType                               flag )
+{
+  for( uint_t k = 0; k < function.getNumberOfBlocks(); k++ ) {
+    createVectorFromFunction( function[k], numerator[k], vec, level, flag ); 
+  }
+}
+
+inline void createFunctionFromVector( const BlockFunction< PetscReal >&   function,
+                                      const BlockFunction< PetscInt >&    numerator,
+                                      const std::shared_ptr< VectorProxy >& vec,
+                                      uint_t                                level,
+                                      DoFType                               flag )
+{
+  for( uint_t k = 0; k < function.getNumberOfBlocks(); k++ ) {
+    createFunctionFromVector( function[k], numerator[k], vec, level, flag ); 
+  }
+}
+
+} // namespace petsc
+} // namespace hyteg
+#endif

src/hyteg/functions/FunctionWrapper.hpp

@@ -24,6 +24,18 @@
 
 #include "hyteg/functions/FunctionTraits.hpp"
 #include "hyteg/functions/GenericFunction.hpp"
+#include "hyteg/sparseassembly/VectorProxy.hpp"
+
+// A whole lot of includes, so that createVectorFromFunction below has
+// a valid prototype for all possible cases
+#include "hyteg/dgfunctionspace/DGPetsc.hpp"
+#include "hyteg/edgedofspace/EdgeDoFPetsc.hpp"
+#include "hyteg/facedofspace/FaceDoFPetsc.hpp"
+#include "hyteg/p1functionspace/P1Petsc.hpp"
+#include "hyteg/p2functionspace/P2Petsc.hpp"
+
+// only needed for using PetscInt in to/fromVector() below!
+#include "hyteg/petsc/PETScWrapper.hpp"
 
 namespace hyteg {
 
@@ -172,6 +184,43 @@ class FunctionWrapper final : public GenericFunction< typename FunctionTrait< fu
       return numberOfLocalDoFs< typename FunctionTrait< WrappedFuncType >::Tag >( *storage, level );
    }
 
+#ifdef HYTEG_BUILD_WITH_PETSC
+   /// conversion to/from linear algebra representation
+   /// @{
+   void toVector( const GenericFunction< PetscInt >&    numerator,
+                  const std::shared_ptr< VectorProxy >& vec,
+                  uint_t                                level,
+                  DoFType                               flag ) const
+   {
+      if constexpr ( std::is_same< value_t, PetscReal >::value )
+      {
+         using numer_t = typename func_t::template FunctionType< PetscInt >;
+         petsc::createVectorFromFunction( *wrappedFunc_, numerator.template unwrap< numer_t >(), vec, level, flag );
+      }
+      else
+      {
+         WALBERLA_ABORT( "FunctionWrapper::toVector() only works for ValueType being identical to PetscReal" );
+      }
+   };
+
+   void fromVector( const GenericFunction< PetscInt >&    numerator,
+                    const std::shared_ptr< VectorProxy >& vec,
+                    uint_t                                level,
+                    DoFType                               flag ) const
+   {
+      if constexpr ( std::is_same< value_t, PetscReal >::value )
+      {
+         using numer_t = typename func_t::template FunctionType< PetscInt >;
+         petsc::createFunctionFromVector( *wrappedFunc_, numerator.template unwrap< numer_t >(), vec, level, flag );
+      }
+      else
+      {
+         WALBERLA_ABORT( "FunctionWrapper::fromVector() only works for ValueType being identical to PetscReal" );
+      }
+   };
+      /// @}
+#endif
+
  private:
    std::unique_ptr< func_t > wrappedFunc_;
 };

src/hyteg/functions/GenericFunction.hpp

@@ -23,7 +23,9 @@
 #include "core/DataTypes.h"
 
 #include "hyteg/boundary/BoundaryConditions.hpp"
+#include "hyteg/petsc/PETScWrapper.hpp"
 #include "hyteg/primitivestorage/PrimitiveStorage.hpp"
+#include "hyteg/sparseassembly/VectorProxy.hpp"
 #include "hyteg/types/flags.hpp"
 
 namespace hyteg {
@@ -125,6 +127,21 @@ class GenericFunction
    virtual void enumerate( uint_t level ) const = 0;
 
    virtual void enumerate( uint_t level, value_t& offset ) const = 0;
+
+#ifdef HYTEG_BUILD_WITH_PETSC
+   /// conversion to/from linear algebra representation
+   /// @{
+   virtual void toVector( const GenericFunction< PetscInt >&    numerator,
+                          const std::shared_ptr< VectorProxy >& vec,
+                          uint_t                                level,
+                          DoFType                               flag ) const = 0;
+
+   virtual void fromVector( const GenericFunction< PetscInt >&    numerator,
+                            const std::shared_ptr< VectorProxy >& vec,
+                            uint_t                                level,
+                            DoFType                               flag ) const = 0;
+   /// @}
+#endif
 };
 
 // Special version of numberOfLocalDoFs for GenericFunctions

src/hyteg/functions/GenericFunctionPetsc.hpp

+/*
+ * Copyright (c) 2021 Marcus Mohr.
+ *
+ * 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/>.
+ */
+
+#pragma once
+
+#include "core/DataTypes.h"
+
+#include "hyteg/functions/GenericFunction.hpp"
+#include "hyteg/petsc/PETScWrapper.hpp"
+#include "hyteg/sparseassembly/VectorProxy.hpp"
+#include "hyteg/types/flags.hpp"
+
+#ifdef HYTEG_BUILD_WITH_PETSC
+namespace hyteg {
+namespace petsc {
+
+inline void createVectorFromFunction( const GenericFunction< PetscReal >&   function,
+                                      const GenericFunction< PetscInt >&    numerator,
+                                      const std::shared_ptr< VectorProxy >& vec,
+                                      uint_t                                level,
+                                      DoFType                               flag )
+{
+   function.toVector( numerator, vec, level, flag );
+}
+
+inline void createFunctionFromVector( const GenericFunction< PetscReal >&   function,
+                                      const GenericFunction< PetscInt >&    numerator,
+                                      const std::shared_ptr< VectorProxy >& vec,
+                                      uint_t                                level,
+                                      DoFType                               flag )
+{
+   function.fromVector( numerator, vec, level, flag );
+}
+
+} // namespace petsc
+} // namespace hyteg
+#endif

src/hyteg/petsc/PETScVector.hpp

@@ -114,6 +114,7 @@ class PETScVector
  protected:
    MPI_Comm petscCommunicator_;
    Vec      vec;
+
 };
 
 } // namespace hyteg

tests/hyteg/BlockFunctions/BlockFunctionBasicTest.cpp

@@ -30,7 +30,10 @@
 #include "hyteg/composites/P2P1TaylorHoodStokesOperator.hpp"
 #include "hyteg/dataexport/VTKOutput.hpp"
 #include "hyteg/functions/BlockFunction.hpp"
+#include "hyteg/functions/BlockFunctionPetsc.hpp"
 #include "hyteg/p1functionspace/P1Function.hpp"
+#include "hyteg/petsc/PETScManager.hpp"
+#include "hyteg/petsc/PETScVector.hpp"
 #include "hyteg/primitivestorage/SetupPrimitiveStorage.hpp"
 #include "hyteg/primitivestorage/loadbalancing/SimpleBalancer.hpp"
 
@@ -191,6 +194,51 @@ void testEnumerate()
    WALBERLA_CHECK_EQUAL( nDoFs, value + 1 );
 }
 
+#ifdef HYTEG_BUILD_WITH_PETSC
+template < template < typename > class FunctionKind >
+void testPETScConversion( const std::string& kind, bool exportVTU = false )
+{
+   WALBERLA_LOG_INFO_ON_ROOT( "RUNNING WITH '" << kind << "'" );
+   // WALBERLA_LOG_INFO_ON_ROOT( "-> Running test for " << FunctionTrait< FunctionKind< real_t > >::getTypeName() );
+
+   MeshInfo              mesh = MeshInfo::fromGmshFile( "../../data/meshes/tri_1el.msh" );
+   SetupPrimitiveStorage setupStorage( mesh, uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) );
+   setupStorage.setMeshBoundaryFlagsOnBoundary( 1, 0, true );
+   std::shared_ptr< walberla::WcTimingTree > timingTree( new walberla::WcTimingTree() );
+   std::shared_ptr< PrimitiveStorage >       storage = std::make_shared< PrimitiveStorage >( setupStorage, timingTree );
+
+   uint_t level = 4;
+
+   FunctionKind< real_t >   src( "src", storage, level, level );
+   FunctionKind< real_t >   dst( "dst", storage, level, level );
+   FunctionKind< PetscInt > numerator( "numerator", storage, level, level );
+
+   std::function< real_t( const hyteg::Point3D& ) > expression = []( const hyteg::Point3D& x ) {
+      real_t value;
+      value = std::sin( real_c( 3 ) * x[0] ) + real_c( 0.5 ) * x[1] * x[1];
+      return value;
+   };
+
+   numerator.enumerate( level );
+   src.interpolate( expression, level );
+
+   PETScVector< real_t, BlockFunction > vector( src, numerator, level, All );
+   vector.createFunctionFromVector( dst, numerator, level, All );
+
+   if ( exportVTU )
+   {
+      std::string fPath = "../../output";
+      std::string fName = "BlockFunctionConversion";
+      WALBERLA_LOG_INFO_ON_ROOT( "Exporting to '" << fPath << "/" << fName << "'" );
+      VTKOutput vtkOutput( fPath, fName, storage );
+      vtkOutput.add( src );
+      vtkOutput.add( dst );
+      vtkOutput.add( numerator );
+      vtkOutput.write( level );
+   }
+}
+#endif
+
 int main( int argc, char* argv[] )
 {
    walberla::debug::enterTestMode();
@@ -209,5 +257,15 @@ int main( int argc, char* argv[] )
    // test enumerating using a P2P1TaylorHoodBlockFunction
    testEnumerate();
 
+#ifdef HYTEG_BUILD_WITH_PETSC
+   WALBERLA_LOG_INFO_ON_ROOT( "" );
+   WALBERLA_LOG_INFO_ON_ROOT( "--------------------------------" );
+   WALBERLA_LOG_INFO_ON_ROOT( " Function <-> Vector Conversion" );
+   WALBERLA_LOG_INFO_ON_ROOT( "--------------------------------" );
+   PETScManager petscManager( &argc, &argv );
+   testPETScConversion< P1P1TaylorHoodStokesBlockFunction >( "P1P1TaylorHoodStokesBlockFunction" );
+   testPETScConversion< P2P1TaylorHoodBlockFunction >( "P2P1TaylorHoodBlockFunction" );
+#endif
+
    return EXIT_SUCCESS;
 }

tests/hyteg/BlockFunctions/FunctionWrapperTest.cpp

@@ -18,6 +18,8 @@
  * along with this program. If not, see <http://www.gnu.org/licenses/>.
  */
 
+#include "hyteg/functions/FunctionWrapper.hpp"
+
 #include "core/DataTypes.h"
 #include "core/Environment.h"
 #include "core/debug/CheckFunctions.h"
@@ -26,32 +28,59 @@
 
 #include "hyteg/facedofspace/FaceDoFFunction.hpp"
 #include "hyteg/functions/FunctionProperties.hpp"
+#include "hyteg/functions/GenericFunctionPetsc.hpp"
 #include "hyteg/mesh/MeshInfo.hpp"
 #include "hyteg/p1functionspace/P1Function.hpp"
 #include "hyteg/p1functionspace/P1VectorFunction.hpp"
 #include "hyteg/p2functionspace/P2Function.hpp"
 #include "hyteg/p2functionspace/P2VectorFunction.hpp"
+#include "hyteg/petsc/PETScManager.hpp"
+#include "hyteg/petsc/PETScVector.hpp"
 #include "hyteg/primitivestorage/SetupPrimitiveStorage.hpp"
 
 // Perform some basic test on the FunctionWrapper class
 
-#include "hyteg/functions/FunctionWrapper.hpp"
-
 using namespace hyteg;
 
 template < typename func_t >
-void wrapFunction( const std::shared_ptr< PrimitiveStorage >& storage,
-                   size_t                                     minLevel,
-                   size_t                                     maxLevel,
-                   uint_t dim )
+void wrapFunction( const std::shared_ptr< PrimitiveStorage >& storage, size_t minLevel, size_t maxLevel, uint_t dim )
 
 {
-  std::string kind = FunctionTrait< func_t >::getTypeName();
-  FunctionWrapper< func_t > wrappedFunc( kind, storage, minLevel, maxLevel );
-  WALBERLA_LOG_INFO_ON_ROOT( "-> wrapped a '" << kind << "'" );
-  WALBERLA_CHECK_EQUAL( dim, wrappedFunc.getDimension() );
+   std::string               kind = FunctionTrait< func_t >::getTypeName();
+   FunctionWrapper< func_t > wrappedFunc( kind, storage, minLevel, maxLevel );
+   WALBERLA_LOG_INFO_ON_ROOT( "-> wrapped a '" << kind << "'" );
+   WALBERLA_CHECK_EQUAL( dim, wrappedFunc.getDimension() );
 }
- 
+
+#ifdef HYTEG_BUILD_WITH_PETSC
+template < template < typename > class FunctionKind >
+void testPETScConversion( const std::shared_ptr< PrimitiveStorage >& storage )
+{
+   WALBERLA_LOG_INFO_ON_ROOT( "-> Running test for " << FunctionTrait< FunctionKind< real_t > >::getTypeName() );
+
+   uint_t level = 4;
+
+   FunctionWrapper< FunctionKind< real_t > >   src( "testing", storage, level, level );
+   FunctionWrapper< FunctionKind< real_t > >   dst( "testing", storage, level, level );
+   FunctionWrapper< FunctionKind< PetscInt > > numerator( "numerator", storage, level, level );
+
+   std::function< real_t( const hyteg::Point3D& ) > expression = []( const hyteg::Point3D& x ) {
+      real_t value;
+      value = std::sin( real_c( 3 ) * x[0] ) + real_c( 0.5 ) * x[1] * x[1];
+      return value;
+   };
+
+   numerator.enumerate( level );
+   src.interpolate( expression, level );
+
+   PETScVector< real_t, GenericFunction > vector( src, numerator, level, All );
+   vector.createFunctionFromVector( dst, numerator, level, All );
+   dst.assign( {real_c( 1 ), real_c( -1 )}, {dst, src}, level, All );
+   // real_t diff = dst.getMaxComponentMagnitude( level, All );
+   // WALBERLA_CHECK_FLOAT_EQUAL( diff, real_c( 0 ) );
+   // WALBERLA_LOG_INFO_ON_ROOT( "   max. difference = " << diff );
+}
+#endif
 
 int main( int argc, char* argv[] )
 {
@@ -171,5 +200,17 @@ int main( int argc, char* argv[] )
    WALBERLA_CHECK_FLOAT_EQUAL( p2vec[0].getMaxMagnitude( maxLevel ), real_c( 4 ) );
    WALBERLA_LOG_INFO_ON_ROOT( "P2VecFunc.interpolate() -> check" );
 
+#ifdef HYTEG_BUILD_WITH_PETSC
+   WALBERLA_LOG_INFO_ON_ROOT( "--------------------------------" );
+   WALBERLA_LOG_INFO_ON_ROOT( " Function <-> Vector Conversion" );
+   WALBERLA_LOG_INFO_ON_ROOT( "--------------------------------" );
+   PETScManager petscManager( &argc, &argv );
+   testPETScConversion< P1Function >( storage );
+   testPETScConversion< P2Function >( storage );
+   testPETScConversion< EdgeDoFFunction >( storage );
+   testPETScConversion< P1VectorFunction >( storage );
+   testPETScConversion< P2VectorFunction >( storage );
+#endif
+
    return EXIT_SUCCESS;
 }

tests/hyteg/forms/SingleRowIntegrationTest.cpp

@@ -70,13 +70,13 @@ real_t
 template < typename Form, uint_t dim, uint_t rows, uint_t cols >
 void compareRows( const Form& form, const std::array< Point3D, dim + 1 >& element, uint_t row, real_t tol )
 {
-   typedef Matrixr< rows, cols > MatType;
-   typedef Matrixr< 1, cols >    RowType;
+   typedef Matrixr< rows, cols > LocalMatType;
+   typedef Matrixr< 1, cols >    LocalRowType;
 
-   MatType elMat;
-   RowType elRow;
-   RowType elRowOfFullMat;
-   RowType elRowDifference;
+   LocalMatType elMat;
+   LocalRowType elRow;
+   LocalRowType elRowOfFullMat;
+   LocalRowType elRowDifference;
 
    form.integrateAll( element, elMat );
    form.integrateRow( row, element, elRow );