diff --git a/src/TNL/Solvers/PDE/ExplicitUpdater.h b/src/TNL/Solvers/PDE/ExplicitUpdater.h
index d6683b79d54a60c39515f3f231526978a60a3852..0559018f36717698a126738629aba970208b7b2c 100644
--- a/src/TNL/Solvers/PDE/ExplicitUpdater.h
+++ b/src/TNL/Solvers/PDE/ExplicitUpdater.h
@@ -120,13 +120,19 @@ class ExplicitUpdater
                                                  typename MeshFunction::DeviceType,
                                                  typename MeshFunction::IndexType > >::value != true,
             "Error: I am getting Vector instead of MeshFunction or similar object. You might forget to bind DofVector into MeshFunction in you method getExplicitUpdate."  );
+         TNL_ASSERT_GT( uPointer->getData().getSize(), 0, "The first MeshFunction in the parameters was not bound." );
+         TNL_ASSERT_GT( fuPointer->getData().getSize(), 0, "The second MeshFunction in the parameters was not bound." );
+
+         TNL_ASSERT_EQ( uPointer->getData().getSize(), mesh.template getEntitiesCount< EntityType >(),
+                        "The first MeshFunction in the parameters was not bound properly." );
+         TNL_ASSERT_EQ( fuPointer->getData().getSize(), mesh.template getEntitiesCount< EntityType >(),
+                        "The second MeshFunction in the parameters was not bound properly." );
             
          TNL_ASSERT_TRUE( this->userData.differentialOperator,
                           "The differential operator is not correctly set-up. Use method setDifferentialOperator() to do it." );
          TNL_ASSERT_TRUE( this->userData.rightHandSide,
                           "The right-hand side is not correctly set-up. Use method setRightHandSide() to do it." );
-         
-         
+
          this->userData.time = time;
          this->userData.u = &uPointer.template modifyData< DeviceType >();
          this->userData.fu = &fuPointer.template modifyData< DeviceType >();
@@ -135,8 +141,6 @@ class ExplicitUpdater
                                                          TraverserInteriorEntitiesProcessor >
                                                        ( meshPointer,
                                                          userData );
-         this->userData.time = time + tau;
-         
       }
       
       template< typename EntityType >
diff --git a/tests/benchmarks/heat-equation-benchmark/HeatEquationBenchmarkProblem_impl.h b/tests/benchmarks/heat-equation-benchmark/HeatEquationBenchmarkProblem_impl.h
index 76539e22f7cce06f985c1406aea3c7aabf356f2e..00fdbf68e7d437753ba0d56054c39f1068d5736f 100644
--- a/tests/benchmarks/heat-equation-benchmark/HeatEquationBenchmarkProblem_impl.h
+++ b/tests/benchmarks/heat-equation-benchmark/HeatEquationBenchmarkProblem_impl.h
@@ -11,6 +11,11 @@
 #include "Tuning/SimpleCell.h"
 #include "Tuning/GridTraverser.h"
 
+//#define WITH_TNL  // In the 'tunning' part, this serves for comparison of performance 
+                  // when using common TNL structures compared to the benchmark ones
+
+
+
 template< typename Mesh,
           typename BoundaryCondition,
           typename RightHandSide,
@@ -534,19 +539,28 @@ getExplicitUpdate( const RealType& time,
       {
          if( std::is_same< DeviceType, Devices::Cuda >::value )
          {   
-            /*using ExplicitUpdaterType = TNL::Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, DifferentialOperator, BoundaryCondition, RightHandSide >;
+            this->u->bind( mesh, uDofs );
+            this->fu->bind( mesh, fuDofs );                     
+            
+            
+            /*this->explicitUpdater.template update< typename Mesh::Cell >( time, tau, mesh, this->u, this->fu );
+            return;*/
+            
+#ifdef WITH_TNL
+            using ExplicitUpdaterType = TNL::Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, DifferentialOperator, BoundaryCondition, RightHandSide >;
             using Cell = typename MeshType::Cell;
             using MeshTraverserType = Meshes::Traverser< MeshType, Cell >;
             using UserData = TNL::Solvers::PDE::ExplicitUpdaterTraverserUserData< RealType,
                MeshFunctionType,
                DifferentialOperator,
                BoundaryCondition,
-               RightHandSide >;*/
+               RightHandSide >;
             
+#else
             //using CellConfig = Meshes::GridEntityNoStencilStorage;
-            
             using CellConfig = Meshes::GridEntityCrossStencilStorage< 1 >;
             using ExplicitUpdaterType = ExplicitUpdater< Mesh, MeshFunctionType, DifferentialOperator, BoundaryCondition, RightHandSide >;
+            //using Cell = typename MeshType::Cell; 
             using Cell = SimpleCell< Mesh, CellConfig >;
             using MeshTraverserType = Traverser< MeshType, Cell >;
             using UserData = ExplicitUpdaterTraverserUserData< RealType,
@@ -554,24 +568,22 @@ getExplicitUpdate( const RealType& time,
                DifferentialOperator,
                BoundaryCondition,
                RightHandSide >;
+#endif            
 
             using InteriorEntitiesProcessor = typename ExplicitUpdaterType::TraverserInteriorEntitiesProcessor;
             using BoundaryEntitiesProcessor = typename ExplicitUpdaterType::TraverserBoundaryEntitiesProcessor;
-         
-            this->u->bind( mesh, uDofs );
-            this->fu->bind( mesh, fuDofs );         
             
             UserData userData;
             userData.time = time;
             userData.differentialOperator = &this->differentialOperatorPointer.template getData< Devices::Cuda >();
             userData.boundaryConditions = &this->boundaryConditionPointer.template getData< Devices::Cuda >();
             userData.rightHandSide = &this->rightHandSidePointer.template getData< Devices::Cuda >();
-            //userData.u = &this->u.template modifyData< Devices::Cuda >(); //uDofs->getData();
-            //userData.fu = &this->fu.template modifyData< Devices::Cuda >(); //fuDofs->getData();
-            userData.u = uDofs->getData();
-            userData.fu = fuDofs->getData();
-            
-                        
+            userData.u = &this->u.template modifyData< Devices::Cuda >(); //uDofs->getData();
+            userData.fu = &this->fu.template modifyData< Devices::Cuda >(); //fuDofs->getData();
+#ifndef WITH_TNL
+            userData.real_u = uDofs->getData();
+            userData.real_fu = fuDofs->getData();
+#endif                        
             const IndexType gridXSize = mesh->getDimensions().x();
             const IndexType gridYSize = mesh->getDimensions().y();
             dim3 cudaBlockSize( 16, 16 );
@@ -580,41 +592,42 @@ getExplicitUpdate( const RealType& time,
             
             TNL::Devices::Cuda::synchronizeDevice();
             int cudaErr;
-            /*Meshes::Traverser< MeshType, Cell > meshTraverser;
-            meshTraverser.template processBoundaryEntities< UserData,
-                                                      BoundaryEntitiesProcessor >
+            Meshes::Traverser< MeshType, Cell > meshTraverser;
+            meshTraverser.template processInteriorEntities< UserData,
+                                                      InteriorEntitiesProcessor >
                                                           ( mesh,
                                                             userData );
-             */
-            
-            
-            _boundaryConditionsKernel< BoundaryEntitiesProcessor, UserData, MeshType, RealType, IndexType >
+             // */
+            /*_heatEquationKernel< InteriorEntitiesProcessor, UserData, MeshType, RealType, IndexType >
             <<< cudaGridSize, cudaBlockSize >>>
                ( &mesh.template getData< Devices::Cuda >(),
                 userData );
-                //&userDataPtr.template modifyData< Devices::Cuda >() );
-            if( ( cudaErr = cudaGetLastError() ) != cudaSuccess )
+                //&userDataPtr.template modifyData< Devices::Cuda >() );*/
+            if( cudaGetLastError() != cudaSuccess )
             {
-               std::cerr << "Setting of boundary conditions failed. " << cudaErr << std::endl;
+               std::cerr << "Laplace operator failed." << std::endl;
                return;
             }
-
             
-            /*meshTraverser.template processInteriorEntities< UserData,
-                                                      InteriorEntitiesProcessor >
+            meshTraverser.template processBoundaryEntities< UserData,
+                                                      BoundaryEntitiesProcessor >
                                                           ( mesh,
                                                             userData );
-             * */
-            _heatEquationKernel< InteriorEntitiesProcessor, UserData, MeshType, RealType, IndexType >
+            // */
+           /*_boundaryConditionsKernel< BoundaryEntitiesProcessor, UserData, MeshType, RealType, IndexType >
             <<< cudaGridSize, cudaBlockSize >>>
                ( &mesh.template getData< Devices::Cuda >(),
                 userData );
-                //&userDataPtr.template modifyData< Devices::Cuda >() );*/
-            if( cudaGetLastError() != cudaSuccess )
+                //&userDataPtr.template modifyData< Devices::Cuda >() );
+            // */ 
+            if( ( cudaErr = cudaGetLastError() ) != cudaSuccess )
             {
-               std::cerr << "Laplace operator failed." << std::endl;
+               std::cerr << "Setting of boundary conditions failed. " << cudaErr << std::endl;
                return;
             }
+
+            
+            
          }
       }      
    }
@@ -629,11 +642,74 @@ HeatEquationBenchmarkProblem< Mesh, BoundaryCondition, RightHandSide, Differenti
 applyBoundaryConditions( const RealType& time,
                          DofVectorPointer& uDofs )
 {
+   const MeshPointer& mesh = this->getMesh();
    if( this->cudaKernelType == "templated" )
    {
       this->bindDofs( uDofs );
       this->explicitUpdater.template applyBoundaryConditions< typename Mesh::Cell >( this->getMesh(), time, this->u );
    }
+   if( this->cudaKernelType == "tunning" )
+   {
+      return;
+      this->bindDofs( uDofs );
+      this->explicitUpdater.template applyBoundaryConditions< typename Mesh::Cell >( this->getMesh(), time, this->u );
+      return;
+      
+#ifdef WITH_TNL
+      using ExplicitUpdaterType = TNL::Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, DifferentialOperator, BoundaryCondition, RightHandSide >;
+      using Cell = typename MeshType::Cell;
+      using MeshTraverserType = Meshes::Traverser< MeshType, Cell >;
+      using UserData = TNL::Solvers::PDE::ExplicitUpdaterTraverserUserData< RealType,
+         MeshFunctionType,
+         DifferentialOperator,
+         BoundaryCondition,
+         RightHandSide >;
+            
+#else
+      //using CellConfig = Meshes::GridEntityNoStencilStorage;
+      using CellConfig = Meshes::GridEntityCrossStencilStorage< 1 >;
+      using ExplicitUpdaterType = ExplicitUpdater< Mesh, MeshFunctionType, DifferentialOperator, BoundaryCondition, RightHandSide >;
+      //using Cell = typename MeshType::Cell; 
+      using Cell = SimpleCell< Mesh, CellConfig >;
+      using MeshTraverserType = Traverser< MeshType, Cell >;
+      using UserData = ExplicitUpdaterTraverserUserData< RealType,
+         MeshFunctionType,
+         DifferentialOperator,
+         BoundaryCondition,
+         RightHandSide >;
+#endif            
+         using InteriorEntitiesProcessor = typename ExplicitUpdaterType::TraverserInteriorEntitiesProcessor;
+         using BoundaryEntitiesProcessor = typename ExplicitUpdaterType::TraverserBoundaryEntitiesProcessor;
+
+         UserData userData;
+         userData.time = time;
+         userData.differentialOperator = &this->differentialOperatorPointer.template getData< Devices::Cuda >();
+         userData.rightHandSide = &this->rightHandSidePointer.template getData< Devices::Cuda >();
+         userData.u = &this->u.template modifyData< Devices::Cuda >(); //uDofs->getData();
+#ifndef WITH_TNL
+         userData.real_u = uDofs->getData();
+#endif
+      userData.boundaryConditions = &this->boundaryConditionPointer.template getData< Devices::Cuda >();
+      Meshes::Traverser< MeshType, Cell > meshTraverser;
+      /*meshTraverser.template processBoundaryEntities< UserData,
+                                                BoundaryEntitiesProcessor >
+                                                    ( mesh,
+                                                      userData );*/
+      // */
+      /*_boundaryConditionsKernel< BoundaryEntitiesProcessor, UserData, MeshType, RealType, IndexType >
+      <<< cudaGridSize, cudaBlockSize >>>
+         ( &mesh.template getData< Devices::Cuda >(),
+          userData );
+          //&userDataPtr.template modifyData< Devices::Cuda >() );
+      // */ 
+      int cudaErr;
+      if( ( cudaErr = cudaGetLastError() ) != cudaSuccess )
+      {
+         std::cerr << "Setting of boundary conditions failed. " << cudaErr << std::endl;
+         return;
+      }
+      
+   }
 }
 
 
diff --git a/tests/benchmarks/heat-equation-benchmark/Tuning/ExplicitUpdater.h b/tests/benchmarks/heat-equation-benchmark/Tuning/ExplicitUpdater.h
index f0d695bf90129ed8f41961d0522028d04d31f2e2..62f0c55c903a83eb5749d73a6e36ea59cdeef8a9 100644
--- a/tests/benchmarks/heat-equation-benchmark/Tuning/ExplicitUpdater.h
+++ b/tests/benchmarks/heat-equation-benchmark/Tuning/ExplicitUpdater.h
@@ -41,9 +41,9 @@ class ExplicitUpdaterTraverserUserData
 
       const RightHandSide* rightHandSide;
 
-      //MeshFunction *uMf, *fuMf;
+      MeshFunction *u, *fu;
       
-      Real *u, *fu;
+      Real *real_u, *real_fu;
       
       ExplicitUpdaterTraverserUserData()
       : time( 0.0 ),
@@ -117,18 +117,16 @@ class ExplicitUpdater
          
          
          this->userDataPointer->time = time;
-         this->userDataPointer->u = uPointer->getData().getData();
-         this->userDataPointer->fu = fuPointer->getData().getData();
+         this->userDataPointer->u = &uPointer->template modifyData< DeviceType >();
+         this->userDataPointer->fu = &fuPointer->template modifyData< DeviceType >();
+         this->userDataPointer->real_u = uPointer->getData().getData();
+         this->userDataPointer->real_fu = fuPointer->getData().getData();         
          TNL::Traverser< MeshType, EntityType > meshTraverser;
          meshTraverser.template processInteriorEntities< TraverserUserData,
                                                          TraverserInteriorEntitiesProcessor >
                                                        ( meshPointer,
                                                          userDataPointer );
          this->userDataPointer->time = time + tau;
-         /*meshTraverser.template processBoundaryEntities< TraverserUserData,
-                                             TraverserBoundaryEntitiesProcessor >
-                                           ( meshPointer,
-                                             userDataPointer );*/  
       }
       
       template< typename EntityType >
@@ -156,19 +154,18 @@ class ExplicitUpdater
                                               TraverserUserData& userData,
                                               const GridEntity& entity )
             {
-               /*( *userData.u )( entity ) = ( *userData.boundaryConditions )
-                  ( *userData.u, entity, userData.time );*/
+               ( *userData.u )( entity ) = ( *userData.boundaryConditions )
+                  ( *userData.u, entity, userData.time );
             }
             
-            //template< typename EntityType >            
+
             __cuda_callable__
             static inline void processEntity( const MeshType& mesh,
                                               TraverserUserData& userData,
-                                              //const EntityType& entity,
                                               const IndexType& entityIndex,
                                               const typename MeshType::CoordinatesType& coordinates )
             {
-               userData.u[ entityIndex ] = 0.0; /*( *userData.boundaryConditions )
+               userData.real_u[ entityIndex ] = 0.0; /* ( *userData.boundaryConditions )
                   ( *userData.u, entity, userData.time );*/
             }
             
@@ -188,23 +185,21 @@ class ExplicitUpdater
                                               const EntityType& entity )
             {
                typedef Functions::FunctionAdapter< MeshType, RightHandSide > FunctionAdapter;
-               ( userData.fu )[ entity.getIndex() ]  = 
-                  ( *userData.differentialOperator )( userData.u, entity, userData.time );
+               ( *userData.fu )( entity )  = 
+                  ( *userData.differentialOperator )( *userData.u, entity, userData.time );
                    + FunctionAdapter::getValue( *userData.rightHandSide, entity, userData.time );
                
             }
 
-            //template< typename EntityType >
             __cuda_callable__
             static inline void processEntity( const MeshType& mesh,
                                               TraverserUserData& userData,
-                                              //const EntityType& entity,
                                               const IndexType& entityIndex,
                                               const typename MeshType::CoordinatesType& coordinates )
             {
                typedef Functions::FunctionAdapter< MeshType, RightHandSide > FunctionAdapter;
-               userData.fu[ entityIndex ] = 
-                       ( *userData.differentialOperator )( mesh, userData.u, entityIndex, coordinates, userData.time );
+               userData.real_fu[ entityIndex ] = 
+                       ( *userData.differentialOperator )( mesh, userData.real_u, entityIndex, coordinates, userData.time );
                     //   + 0.0;
             }
             
diff --git a/tests/benchmarks/heat-equation-benchmark/Tuning/tunning.h b/tests/benchmarks/heat-equation-benchmark/Tuning/tunning.h
index 4c5384345b64329fbe15118b3cfa91a026d9b083..07a5d7f62bdd7e08fb66a7df0f1e9c8a869cb17b 100644
--- a/tests/benchmarks/heat-equation-benchmark/Tuning/tunning.h
+++ b/tests/benchmarks/heat-equation-benchmark/Tuning/tunning.h
@@ -31,6 +31,9 @@ struct Data
 
 #ifdef HAVE_CUDA
 
+#define WITH_CELL  // Serves for comparison of performance when using SimpleCell
+                   // vs. using only cell index and coordinates
+
 template< typename BoundaryEntitiesProcessor, typename UserData, typename Grid, typename Real, typename Index >
 __global__ void _boundaryConditionsKernel( const Grid* grid,
                                            UserData userData )
@@ -40,30 +43,56 @@ __global__ void _boundaryConditionsKernel( const Grid* grid,
    using Coordinates = typename Grid::CoordinatesType;
    const Index& gridXSize = grid->getDimensions().x();
    const Index& gridYSize = grid->getDimensions().y();
+#ifdef WITH_CELL   
+   SimpleCell< Grid > cell( *grid,
+      Coordinates( ( blockIdx.x ) * blockDim.x + threadIdx.x,
+                   ( blockIdx.y ) * blockDim.y + threadIdx.y ) );
+   Coordinates& coordinates = cell.getCoordinates();
+   cell.refresh();   
+#else   
    Coordinates coordinates( ( blockIdx.x ) * blockDim.x + threadIdx.x,
                              ( blockIdx.y ) * blockDim.y + threadIdx.y );
-   
    Index entityIndex = coordinates.y() * gridXSize + coordinates.x();
+#endif   
+   
    
    if( coordinates.x() == 0 && coordinates.y() < gridYSize )
    {
       //u[ c ] = ( *bc )( *grid, u, c, coordinates, 0 );
+#ifdef WITH_CELL
+      BoundaryEntitiesProcessor::processEntity( *grid, userData, cell );
+#else
       BoundaryEntitiesProcessor::processEntity( *grid, userData, entityIndex, coordinates );
+#endif 
    }
    if( coordinates.x() == gridXSize - 1 && coordinates.y() < gridYSize )
    {
       //u[ c ] = ( *bc )( *grid, u, c, coordinates, 0 );
+
+#ifdef WITH_CELL
+      BoundaryEntitiesProcessor::processEntity( *grid, userData, cell );
+#else
       BoundaryEntitiesProcessor::processEntity( *grid, userData, entityIndex, coordinates );
+#endif      
    }
    if( coordinates.y() == 0 && coordinates.x() < gridXSize )
    {
       //u[ c ] = ( *bc )( *grid, u, c, coordinates, 0 );
+#ifdef WITH_CELL
+      BoundaryEntitiesProcessor::processEntity( *grid, userData, cell );
+#else
       BoundaryEntitiesProcessor::processEntity( *grid, userData, entityIndex, coordinates );
+#endif      
    }
    if( coordinates.y() == gridYSize -1  && coordinates.x() < gridXSize )
    {
       //u[ c ] = ( *bc )( *grid, u, c, coordinates, 0 );
+
+#ifdef WITH_CELL
+      BoundaryEntitiesProcessor::processEntity( *grid, userData, cell );
+#else
       BoundaryEntitiesProcessor::processEntity( *grid, userData, entityIndex, coordinates );
+#endif      
    }         
 }
 
@@ -93,12 +122,14 @@ __global__ void _heatEquationKernel( const Grid* grid,
    //if( coordinates.x() > 0 && coordinates.x() < gridXSize - 1 &&
    //    coordinates.y() > 0 && coordinates.y() < gridYSize - 1 )      
    {
+#ifdef WITH_CELL      
       cell.refresh();
       InteriorEntitiesProcessor::processEntity( *grid, userData, cell );
-      
+#else      
       //const Index entityIndex = cell.getCoordinates().y() * gridXSize + cell.getCoordinates().x();
-      //const Index entityIndex = coordinates.y() * gridXSize + coordinates.x();
-      //InteriorEntitiesProcessor::processEntity( *grid, userData, entityIndex, coordinates );
+      const Index entityIndex = coordinates.y() * gridXSize + coordinates.x();
+      InteriorEntitiesProcessor::processEntity( *grid, userData, entityIndex, cell.getCoordinates() );
+#endif      
       
       
       //fu[ entityIndex ] = ( *op )( *grid, userData.u, entityIndex, coordinates, userData.time ); // + 0.1;