diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 41198f09bb19f0a74f0f37a907c0322b1533f1d8..f4bc90294f926270f3485ff67c65bdc0affdb789 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -14,3 +14,4 @@ add_subdirectory( inviscid-flow )
 add_subdirectory( inviscid-flow-sw )
 #add_subdirectory( mean-curvature-flow )
 add_subdirectory( flow )
+add_subdirectory( flow-sw )
diff --git a/examples/inviscid-flow-sw/Upwind.h b/examples/inviscid-flow-sw/Upwind.h
index 38776d80a08b0ed9042b241ebef83abc4221becc..263da044a2edaca855b6c6f3fd050bd10cc7c689 100644
--- a/examples/inviscid-flow-sw/Upwind.h
+++ b/examples/inviscid-flow-sw/Upwind.h
@@ -94,6 +94,14 @@ class Upwind
          this->momentumZOperatorPointer->setPressure( pressure );
          this->energyOperatorPointer->setPressure( pressure );
       }
+
+      void setDensity( const MeshFunctionPointer& density )
+      {
+         this->momentumXOperatorPointer->setDensity( density );
+         this->momentumYOperatorPointer->setDensity( density );
+         this->momentumZOperatorPointer->setDensity( density );
+         this->energyOperatorPointer->setDensity( density );
+      }
       
       void setVelocity( const VectorFieldPointer& velocity )
       {
diff --git a/examples/inviscid-flow-sw/UpwindContinuity.h b/examples/inviscid-flow-sw/UpwindContinuity.h
index 19150164f063d8df12f94a4ea9b26ed591bf3487..cc9b84aba02a102a745195e179bffe28c216d7bd 100644
--- a/examples/inviscid-flow-sw/UpwindContinuity.h
+++ b/examples/inviscid-flow-sw/UpwindContinuity.h
@@ -343,23 +343,23 @@ class UpwindContinuity< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, In
          const RealType& velocity_z_up     = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ up ];
          const RealType& velocity_z_down   = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ down ];
          
-         return -hxInverse * ( 
-                                  u[ center ] / ( 2 * this->gamma ) * ( ( 2 * this->gamma - 1 ) * velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                - u[ west ]   / ( 2 * this->gamma ) * ( ( 2 * this->gamma - 1 ) * velocity_x_west   + std::sqrt( this->gamma * pressure_west   / u[ west ]   ) )
-                                - u[ center ] / ( 2 * this->gamma ) * ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                + u[ east ]   / ( 2 * this->gamma ) * ( velocity_x_east   - std::sqrt( this->gamma * pressure_east   / u[ east   ] ) )
+         return -hxInverse * (
+                                   this->positiveDensityFlux( u[ center ], velocity_x_center, pressure_center )
+                                -  this->positiveDensityFlux( u[ west   ], velocity_x_west  , pressure_west   )
+                                -  this->negativeDensityFlux( u[ center ], velocity_x_center, pressure_center )
+                                +  this->negativeDensityFlux( u[ east   ], velocity_x_east  , pressure_east   )
                              )
-                -hyInverse * ( 
-                                  u[ center ] / ( 2 * this->gamma ) * ( ( 2 * this->gamma - 1 ) * velocity_y_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                - u[ south ]  / ( 2 * this->gamma ) * ( ( 2 * this->gamma - 1 ) * velocity_y_south  + std::sqrt( this->gamma * pressure_south  / u[ south ]  ) )
-                                - u[ center ] / ( 2 * this->gamma ) * ( velocity_y_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                + u[ north ]  / ( 2 * this->gamma ) * ( velocity_y_north  - std::sqrt( this->gamma * pressure_north  / u[ north ]  ) )
+                -hyInverse * (
+                                   this->positiveDensityFlux( u[ center ], velocity_y_center, pressure_center )
+                                -  this->positiveDensityFlux( u[ south  ], velocity_y_south , pressure_south  )
+                                -  this->negativeDensityFlux( u[ center ], velocity_y_center, pressure_center )
+                                +  this->negativeDensityFlux( u[ north  ], velocity_y_north , pressure_north  )
                              )
-                -hzInverse * ( 
-                                  u[ center ] / ( 2 * this->gamma ) * ( ( 2 * this->gamma - 1 ) * velocity_z_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                - u[ down ]   / ( 2 * this->gamma ) * ( ( 2 * this->gamma - 1 ) * velocity_z_down   + std::sqrt( this->gamma * pressure_down   / u[ down ]   ) )
-                                - u[ center ] / ( 2 * this->gamma ) * ( velocity_z_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) ) 
-                                + u[ up ]     / ( 2 * this->gamma ) * ( velocity_z_up     - std::sqrt( this->gamma * pressure_up     / u[ up ]     ) ) 
+                -hzInverse * (
+                                   this->positiveDensityFlux( u[ center ], velocity_z_center, pressure_center )
+                                -  this->positiveDensityFlux( u[ down   ], velocity_z_down  , pressure_down   )
+                                -  this->negativeDensityFlux( u[ center ], velocity_z_center, pressure_center )
+                                +  this->negativeDensityFlux( u[ up     ], velocity_z_up    , pressure_up     )
                              );
          
       }
diff --git a/examples/inviscid-flow-sw/UpwindEnergy.h b/examples/inviscid-flow-sw/UpwindEnergy.h
index 44792f8b3c7ae8094b4549a47c939098877dd5d7..246b5d18b070cb26b390b6a46aa946ac3af05d73 100644
--- a/examples/inviscid-flow-sw/UpwindEnergy.h
+++ b/examples/inviscid-flow-sw/UpwindEnergy.h
@@ -63,60 +63,17 @@ class UpwindEnergyBase
       {
           this->pressure = pressure;
       };
+
+      void setDensity( const MeshFunctionPointer& density )
+      {
+          this->density = density;
+      };
       
       void setArtificialViscosity( const RealType& artificialViscosity )
       {
          this->artificialViscosity = artificialViscosity;
       };
 
-      RealType positiveEnergyFlux( const RealType& density, const RealType& velocity_main, const RealType& velocity_other1, const RealType& velocity_other2, const RealType& pressure ) const
-      {
-         const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
-         const RealType& machNumber = velocity_main / speedOfSound;
-         if ( machNumber <= -1.0 )
-            return 0.0;
-        else if ( machNumber <= 0.0 )
-            return density * speedOfSound * speedOfSound * speedOfSound / ( 2.0 * this->gamma ) 
-                 * ( 
-                     ( ( machNumber + 1.0 ) * 0.5 *  ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) )
-                   + ( machNumber * ( machNumber + 1.0 ) )
-                   + ( ( machNumber + 1.0 ) / ( this->gamma - 1.0 ) )
-                   );
-        else if ( machNumber <= 1.0 )
-            return density * speedOfSound * speedOfSound * speedOfSound / ( 2.0 * this->gamma ) 
-                 * ( 
-                     ( ( ( 2.0 * this->gamma - 1.0 ) * machNumber + 1.0 ) * 0.5 * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) )
-                   + ( machNumber * ( machNumber + 1.0 ) )
-                   + ( ( machNumber + 1.0 ) / ( this->gamma - 1.0 ) )
-                   );
-        else   
-            return velocity_main * ( pressure + pressure / ( this->gamma - 1.0 ) + 0.5 * density * ( velocity_main * velocity_main + velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) );
-      };
-
-      RealType negativeEnergyFlux( const RealType& density, const RealType& velocity_main, const RealType& velocity_other1, const RealType& velocity_other2, const RealType& pressure ) const
-      {
-         const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
-         const RealType& machNumber = velocity_main / speedOfSound;
-         if ( machNumber <= -1.0 )
-            return velocity_main * ( pressure + pressure / ( this->gamma - 1.0 ) + 0.5 * density * ( velocity_main * velocity_main + velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) );
-        else if ( machNumber <= 0.0 )
-            return density * speedOfSound * speedOfSound * speedOfSound / ( 2.0 * this->gamma ) 
-                 * ( 
-                     ( ( ( 2.0 * this->gamma - 1.0 ) * machNumber - 1.0 ) * 0.5 * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) )
-                   - ( machNumber * ( machNumber - 1.0 ) )
-                   + ( ( machNumber - 1.0 ) / ( this->gamma - 1.0 ) )
-                   );
-        else if ( machNumber <= 1.0 )
-            return density * speedOfSound * speedOfSound * speedOfSound / ( 2.0 * this->gamma ) 
-                 * ( 
-                     ( ( machNumber - 1.0 ) * 0.5 * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) )
-                   - ( machNumber * ( machNumber - 1.0 ) )
-                   + ( ( machNumber - 1.0 ) / ( this->gamma - 1.0 ) )
-                   );
-        else 
-            return 0.0;
-      };      
-
       protected:
          
          RealType tau;
@@ -128,6 +85,8 @@ class UpwindEnergyBase
          MeshFunctionPointer pressure;
          
          RealType artificialViscosity;
+
+         MeshFunctionPointer density;
 };
    
 template< typename Mesh,
@@ -159,7 +118,55 @@ class UpwindEnergy< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index
       using typename BaseType::VelocityFieldType;
       using typename BaseType::VelocityFieldPointer;
       using BaseType::Dimensions;      
+
+      RealType positiveEnergyFlux( const RealType& density, const RealType& velocity_main, const RealType& pressure ) const
+      {
+         const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
+         const RealType& machNumber = velocity_main / speedOfSound;
+         if ( machNumber <= -1.0 )
+            return 0.0;
+        else if ( machNumber <= 0.0 )
+            return density * speedOfSound * speedOfSound * speedOfSound / ( 2.0 * this->gamma ) 
+                 * ( 
+                     ( ( machNumber + 1.0 ) * 0.5 *  ( machNumber * machNumber / ( speedOfSound * speedOfSound ) ) )
+                   + ( machNumber * ( machNumber + 1.0 ) )
+                   + ( ( machNumber + 1.0 ) / ( this->gamma - 1.0 ) )
+                   );
+        else if ( machNumber <= 1.0 )
+            return density * speedOfSound * speedOfSound * speedOfSound / ( 2.0 * this->gamma ) 
+                 * ( 
+                     ( ( ( 2.0 * this->gamma - 1.0 ) * machNumber + 1.0 ) * 0.5 * ( machNumber * machNumber / ( speedOfSound * speedOfSound ) ) )
+                   + ( machNumber * ( machNumber + 1.0 ) )
+                   + ( ( machNumber + 1.0 ) / ( this->gamma - 1.0 ) )
+                   );
+        else   
+            return velocity_main * ( pressure + pressure / ( this->gamma - 1.0 ) + 0.5 * density * ( velocity_main * velocity_main ) );
+      };
       
+      RealType negativeEnergyFlux( const RealType& density, const RealType& velocity_main, const RealType& pressure ) const
+      {
+         const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
+         const RealType& machNumber = velocity_main / speedOfSound;
+         if ( machNumber <= -1.0 )
+            return velocity_main * ( pressure + pressure / ( this->gamma - 1.0 ) + 0.5 * density * ( velocity_main * velocity_main ) );
+        else if ( machNumber <= 0.0 )
+            return density * speedOfSound * speedOfSound * speedOfSound / ( 2.0 * this->gamma ) 
+                 * ( 
+                     ( ( ( 2.0 * this->gamma - 1.0 ) * machNumber - 1.0 ) * 0.5 * ( machNumber * machNumber / ( speedOfSound * speedOfSound ) ) )
+                   - ( machNumber * ( machNumber - 1.0 ) )
+                   + ( ( machNumber - 1.0 ) / ( this->gamma - 1.0 ) )
+                   );
+        else if ( machNumber <= 1.0 )
+            return density * speedOfSound * speedOfSound * speedOfSound / ( 2.0 * this->gamma ) 
+                 * ( 
+                     ( ( machNumber - 1.0 ) * 0.5 * ( machNumber * machNumber / ( speedOfSound * speedOfSound ) ) )
+                   - ( machNumber * ( machNumber - 1.0 ) )
+                   + ( ( machNumber - 1.0 ) / ( this->gamma - 1.0 ) )
+                   );
+        else 
+            return 0.0;
+      };      
+
       template< typename MeshFunction, typename MeshEntity >
       __cuda_callable__
       Real operator()( const MeshFunction& u,
@@ -180,15 +187,19 @@ class UpwindEnergy< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index
          const RealType& pressure_west   = this->pressure.template getData< DeviceType >()[ west ];
          const RealType& pressure_east   = this->pressure.template getData< DeviceType >()[ east ];
 
+         const RealType& density_center = this->density.template getData< DeviceType >()[ center ];
+         const RealType& density_west   = this->density.template getData< DeviceType >()[ west ];
+         const RealType& density_east   = this->density.template getData< DeviceType >()[ east ];
+
          const RealType& velocity_x_center = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ center ];
          const RealType& velocity_x_east   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ];
          const RealType& velocity_x_west   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
 
          return -hxInverse * ( 
-                                   this->positiveEnergyFlux( u[ center ], velocity_x_center, 0, 0, pressure_center)
-                                 - this->positiveEnergyFlux( u[ west   ], velocity_x_west  , 0, 0, pressure_west  )
-                                 - this->negativeEnergyFlux( u[ center ], velocity_x_center, 0, 0, pressure_center)
-                                 + this->negativeEnergyFlux( u[ east   ], velocity_x_east  , 0, 0, pressure_east  ) 
+                                   this->positiveEnergyFlux( density_center, velocity_x_center, pressure_center)
+                                 - this->positiveEnergyFlux( density_west  , velocity_x_west  , pressure_west  )
+                                 - this->negativeEnergyFlux( density_center, velocity_x_center, pressure_center)
+                                 + this->negativeEnergyFlux( density_east  , velocity_x_east  , pressure_east  ) 
                              );  
   
       }
@@ -232,7 +243,54 @@ class UpwindEnergy< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index
       using typename BaseType::VelocityFieldType;
       using typename BaseType::VelocityFieldPointer;
       using BaseType::Dimensions;
-      
+
+      RealType positiveEnergyFlux( const RealType& density, const RealType& velocity_main, const RealType& velocity_other1, const RealType& pressure ) const
+      {
+         const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
+         const RealType& machNumber = velocity_main / speedOfSound;
+         if ( machNumber <= -1.0 )
+            return 0.0;
+        else if ( machNumber <= 0.0 )
+            return density * speedOfSound * speedOfSound * speedOfSound / ( 2.0 * this->gamma ) 
+                 * ( 
+                     ( ( machNumber + 1.0 ) * 0.5 *  ( machNumber * machNumber + ( velocity_other1 * velocity_other1 ) / ( speedOfSound * speedOfSound ) ) )
+                   + ( machNumber * ( machNumber + 1.0 ) )
+                   + ( ( machNumber + 1.0 ) / ( this->gamma - 1.0 ) )
+                   );
+        else if ( machNumber <= 1.0 )
+            return density * speedOfSound * speedOfSound * speedOfSound / ( 2.0 * this->gamma ) 
+                 * ( 
+                     ( ( ( 2.0 * this->gamma - 1.0 ) * machNumber + 1.0 ) * 0.5 * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 ) / ( speedOfSound * speedOfSound ) ) )
+                   + ( machNumber * ( machNumber + 1.0 ) )
+                   + ( ( machNumber + 1.0 ) / ( this->gamma - 1.0 ) )
+                   );
+        else   
+            return velocity_main * ( pressure + pressure / ( this->gamma - 1.0 ) + 0.5 * density * ( velocity_main * velocity_main + velocity_other1 * velocity_other1 ) );
+      };
+
+      RealType negativeEnergyFlux( const RealType& density, const RealType& velocity_main, const RealType& velocity_other1, const RealType& pressure ) const
+      {
+         const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
+         const RealType& machNumber = velocity_main / speedOfSound;
+         if ( machNumber <= -1.0 )
+            return velocity_main * ( pressure + pressure / ( this->gamma - 1.0 ) + 0.5 * density * ( velocity_main * velocity_main + velocity_other1 * velocity_other1 ) );
+        else if ( machNumber <= 0.0 )
+            return density * speedOfSound * speedOfSound * speedOfSound / ( 2.0 * this->gamma ) 
+                 * ( 
+                     ( ( ( 2.0 * this->gamma - 1.0 ) * machNumber - 1.0 ) * 0.5 * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 ) / ( speedOfSound * speedOfSound ) ) )
+                   - ( machNumber * ( machNumber - 1.0 ) )
+                   + ( ( machNumber - 1.0 ) / ( this->gamma - 1.0 ) )
+                   );
+        else if ( machNumber <= 1.0 )
+            return density * speedOfSound * speedOfSound * speedOfSound / ( 2.0 * this->gamma ) 
+                 * ( 
+                     ( ( machNumber - 1.0 ) * 0.5 * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 ) / ( speedOfSound * speedOfSound ) ) )
+                   - ( machNumber * ( machNumber - 1.0 ) )
+                   + ( ( machNumber - 1.0 ) / ( this->gamma - 1.0 ) )
+                   );
+        else 
+            return 0.0;
+      };  
 
       template< typename MeshFunction, typename MeshEntity >
       __cuda_callable__
@@ -259,6 +317,12 @@ class UpwindEnergy< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index
          const RealType& pressure_north  = this->pressure.template getData< DeviceType >()[ north ];
          const RealType& pressure_south  = this->pressure.template getData< DeviceType >()[ south ];
 
+         const RealType& density_center = this->density.template getData< DeviceType >()[ center ];
+         const RealType& density_west   = this->density.template getData< DeviceType >()[ west ];
+         const RealType& density_east   = this->density.template getData< DeviceType >()[ east ];
+         const RealType& density_north  = this->density.template getData< DeviceType >()[ north ];
+         const RealType& density_south  = this->density.template getData< DeviceType >()[ south ];
+
          const RealType& velocity_x_center = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ center ];
          const RealType& velocity_x_east   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ];
          const RealType& velocity_x_west   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
@@ -272,16 +336,16 @@ class UpwindEnergy< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index
          const RealType& velocity_y_south  = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ];         
          
          return -hxInverse * ( 
-                                   this->positiveEnergyFlux( u[ center ], velocity_x_center, velocity_y_center, 0, pressure_center)
-                                 - this->positiveEnergyFlux( u[ west   ], velocity_x_west  , velocity_y_west  , 0, pressure_west  )
-                                 - this->negativeEnergyFlux( u[ center ], velocity_x_center, velocity_y_center, 0, pressure_center)
-                                 + this->negativeEnergyFlux( u[ east   ], velocity_x_east  , velocity_y_east  , 0, pressure_east  ) 
+                                   this->positiveEnergyFlux( density_center, velocity_x_center, velocity_y_center, pressure_center)
+                                 - this->positiveEnergyFlux( density_west  , velocity_x_west  , velocity_y_west  , pressure_west  )
+                                 - this->negativeEnergyFlux( density_center, velocity_x_center, velocity_y_center, pressure_center)
+                                 + this->negativeEnergyFlux( density_east  , velocity_x_east  , velocity_y_east  , pressure_east  ) 
                              ) 
                 -hyInverse * ( 
-                                   this->positiveEnergyFlux( u[ center ], velocity_y_center, velocity_x_center, 0, pressure_center)
-                                 - this->positiveEnergyFlux( u[ south  ], velocity_y_south , velocity_x_south , 0, pressure_south )
-                                 - this->negativeEnergyFlux( u[ center ], velocity_y_center, velocity_x_center, 0, pressure_center)
-                                 + this->negativeEnergyFlux( u[ north  ], velocity_y_north , velocity_x_north , 0, pressure_north ) 
+                                   this->positiveEnergyFlux( density_center, velocity_y_center, velocity_x_center, pressure_center)
+                                 - this->positiveEnergyFlux( density_south , velocity_y_south , velocity_x_south , pressure_south )
+                                 - this->negativeEnergyFlux( density_center, velocity_y_center, velocity_x_center, pressure_center)
+                                 + this->negativeEnergyFlux( density_north , velocity_y_north , velocity_x_north , pressure_north ) 
                              );     
       }
 
@@ -325,6 +389,54 @@ class UpwindEnergy< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index
       using typename BaseType::VelocityFieldPointer;
       using BaseType::Dimensions;      
 
+      RealType positiveEnergyFlux( const RealType& density, const RealType& velocity_main, const RealType& velocity_other1, const RealType& velocity_other2, const RealType& pressure ) const
+      {
+         const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
+         const RealType& machNumber = velocity_main / speedOfSound;
+         if ( machNumber <= -1.0 )
+            return 0.0;
+        else if ( machNumber <= 0.0 )
+            return density * speedOfSound * speedOfSound * speedOfSound / ( 2.0 * this->gamma ) 
+                 * ( 
+                     ( ( machNumber + 1.0 ) * 0.5 *  ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) )
+                   + ( machNumber * ( machNumber + 1.0 ) )
+                   + ( ( machNumber + 1.0 ) / ( this->gamma - 1.0 ) )
+                   );
+        else if ( machNumber <= 1.0 )
+            return density * speedOfSound * speedOfSound * speedOfSound / ( 2.0 * this->gamma ) 
+                 * ( 
+                     ( ( ( 2.0 * this->gamma - 1.0 ) * machNumber + 1.0 ) * 0.5 * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) )
+                   + ( machNumber * ( machNumber + 1.0 ) )
+                   + ( ( machNumber + 1.0 ) / ( this->gamma - 1.0 ) )
+                   );
+        else   
+            return velocity_main * ( pressure + pressure / ( this->gamma - 1.0 ) + 0.5 * density * ( velocity_main * velocity_main + velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) );
+      };
+
+      RealType negativeEnergyFlux( const RealType& density, const RealType& velocity_main, const RealType& velocity_other1, const RealType& velocity_other2, const RealType& pressure ) const
+      {
+         const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
+         const RealType& machNumber = velocity_main / speedOfSound;
+         if ( machNumber <= -1.0 )
+            return velocity_main * ( pressure + pressure / ( this->gamma - 1.0 ) + 0.5 * density * ( velocity_main * velocity_main + velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) );
+        else if ( machNumber <= 0.0 )
+            return density * speedOfSound * speedOfSound * speedOfSound / ( 2.0 * this->gamma ) 
+                 * ( 
+                     ( ( ( 2.0 * this->gamma - 1.0 ) * machNumber - 1.0 ) * 0.5 * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) )
+                   - ( machNumber * ( machNumber - 1.0 ) )
+                   + ( ( machNumber - 1.0 ) / ( this->gamma - 1.0 ) )
+                   );
+        else if ( machNumber <= 1.0 )
+            return density * speedOfSound * speedOfSound * speedOfSound / ( 2.0 * this->gamma ) 
+                 * ( 
+                     ( ( machNumber - 1.0 ) * 0.5 * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) )
+                   - ( machNumber * ( machNumber - 1.0 ) )
+                   + ( ( machNumber - 1.0 ) / ( this->gamma - 1.0 ) )
+                   );
+        else 
+            return 0.0;
+      };      
+
       template< typename MeshFunction, typename MeshEntity >
       __cuda_callable__
       Real operator()( const MeshFunction& u,
@@ -355,6 +467,14 @@ class UpwindEnergy< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index
          const RealType& pressure_up     = this->pressure.template getData< DeviceType >()[ up ];
          const RealType& pressure_down   = this->pressure.template getData< DeviceType >()[ down ];
          
+         const RealType& density_center = this->density.template getData< DeviceType >()[ center ];
+         const RealType& density_west   = this->density.template getData< DeviceType >()[ west ];
+         const RealType& density_east   = this->density.template getData< DeviceType >()[ east ];
+         const RealType& density_north  = this->density.template getData< DeviceType >()[ north ];
+         const RealType& density_south  = this->density.template getData< DeviceType >()[ south ];
+         const RealType& density_up     = this->density.template getData< DeviceType >()[ up ];
+         const RealType& density_down   = this->density.template getData< DeviceType >()[ down ];
+         
          const RealType& velocity_x_center = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ center ];
          const RealType& velocity_x_east   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ];
          const RealType& velocity_x_west   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
@@ -379,133 +499,24 @@ class UpwindEnergy< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index
          const RealType& velocity_z_up     = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ up ];
          const RealType& velocity_z_down   = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ down ];         
          
-         return -hxInverse * (
-                               ( u[ center ] / ( 2 * this->gamma ) )
-                               * ( 
-                                   ( ( 2 * this->gamma - 1 ) * velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                 * ( velocity_x_center * velocity_x_center + velocity_y_center * velocity_y_center + velocity_z_center * velocity_z_center ) / 2
-                                 + velocity_x_center * std::sqrt( this->gamma * pressure_center / u[ center ] )
-                                 * ( velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                 + this->gamma * pressure_center / u[ center ] 
-                                 * ( velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                 / ( this->gamma - 1 )
-                                 )
-                             - ( u[ west ] / ( 2 * this->gamma ) )
-                               * ( 
-                                   ( ( 2 * this->gamma - 1 ) * velocity_x_west + std::sqrt( this->gamma * pressure_west / u[ west ] ) )
-                                 * ( velocity_x_west * velocity_x_west + velocity_y_west * velocity_y_west + velocity_z_west * velocity_z_west ) / 2
-                                 + velocity_x_west * std::sqrt( this->gamma * pressure_west / u[ west ] )
-                                 * ( velocity_x_west + std::sqrt( this->gamma * pressure_west / u[ west ] ) )
-                                 + this->gamma * pressure_west / u[ west ]
-                                 * ( velocity_x_west + std::sqrt( this->gamma * pressure_west / u[ west ] ) )
-                                 / ( this->gamma - 1 )
-                                 )
-                             - ( u[ center ] / ( 2 * this->gamma ) )
-                               * ( 
-                                   ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                 * ( velocity_x_center * velocity_x_center + velocity_y_center * velocity_y_center + velocity_z_center * velocity_z_center ) / 2
-                                 + velocity_x_center * std::sqrt( this->gamma * pressure_center / u[ center ] )
-                                 * ( std::sqrt( this->gamma * pressure_center / u[ center ] ) - velocity_x_center )
-                                 + this->gamma * pressure_center / u[ center ]
-                                 * ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                 / ( this->gamma - 1 )
-                                 )
-                             + ( u[ east  ] / ( 2 * this->gamma ) )
-                               * ( 
-                                   ( velocity_x_east  - std::sqrt( this->gamma * pressure_east  / u[ east  ] ) )
-                                 * ( velocity_x_east  * velocity_x_east  + velocity_y_east  * velocity_y_east + velocity_z_east * velocity_z_east  ) / 2
-                                 + velocity_x_east  * std::sqrt( this->gamma * pressure_east  / u[ east  ] )
-                                 * ( std::sqrt( this->gamma * pressure_east  / u[ east  ] ) - velocity_x_east  )
-                                 + this->gamma * pressure_east  / u[ east  ]
-                                 * ( velocity_x_east  - std::sqrt( this->gamma * pressure_east  / u[ east  ] ) )
-                                 / ( this->gamma - 1 )
-                                 )
+         return -hxInverse * ( 
+                                   this->positiveEnergyFlux( density_center, velocity_x_center, velocity_y_center, velocity_z_center, pressure_center)
+                                 - this->positiveEnergyFlux( density_west  , velocity_x_west  , velocity_y_west  , velocity_z_west  , pressure_west  )
+                                 - this->negativeEnergyFlux( density_center, velocity_x_center, velocity_y_center, velocity_z_center, pressure_center)
+                                 + this->negativeEnergyFlux( density_east  , velocity_x_east  , velocity_y_east  , velocity_z_east  , pressure_east  ) 
                              ) 
-                -hyInverse * (
-                               ( u[ center ] / ( 2 * this->gamma ) )
-                               * ( 
-                                   ( ( 2 * this->gamma - 1 ) * velocity_y_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                 * ( velocity_x_center * velocity_x_center + velocity_y_center * velocity_y_center + velocity_z_center * velocity_z_center ) / 2
-                                 + velocity_y_center * std::sqrt( this->gamma * pressure_center / u[ center ] )
-                                 * ( velocity_y_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                 + this->gamma * pressure_center / u[ center ] 
-                                 * ( velocity_y_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                 / ( this->gamma - 1 )
-                                 )
-                             - ( u[ south ] / ( 2 * this->gamma ) )
-                               * ( 
-                                   ( ( 2 * this->gamma - 1 ) * velocity_y_south + std::sqrt( this->gamma * pressure_south / u[ south ] ) )
-                                 * ( velocity_x_south * velocity_x_south + velocity_y_south * velocity_y_south + velocity_z_south * velocity_z_south ) / 2
-                                 + velocity_y_south * std::sqrt( this->gamma * pressure_south / u[ south ] )
-                                 * ( velocity_y_south + std::sqrt( this->gamma * pressure_south / u[ south ] ) )
-                                 + this->gamma * pressure_south / u[ south ]
-                                 * ( velocity_y_south + std::sqrt( this->gamma * pressure_south / u[ south ] ) )
-                                 / ( this->gamma - 1 )
-                                 )
-                             - ( u[ center ] / ( 2 * this->gamma ) )
-                               * ( 
-                                   ( velocity_y_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                 * ( velocity_x_center * velocity_x_center + velocity_y_center * velocity_y_center + velocity_z_center * velocity_z_center ) / 2
-                                 + velocity_y_center * std::sqrt( this->gamma * pressure_center / u[ center ] )
-                                 * ( std::sqrt( this->gamma * pressure_center / u[ center ] ) - velocity_y_center )
-                                 + this->gamma * pressure_center / u[ center ]
-                                 * ( velocity_y_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                 / ( this->gamma - 1 )
-                                 )
-                             + ( u[ north ] / ( 2 * this->gamma ) )
-                               * ( 
-                                   ( velocity_y_north - std::sqrt( this->gamma * pressure_north / u[ north ] ) )
-                                 * ( velocity_x_north * velocity_x_north + velocity_y_north * velocity_y_north + velocity_z_north * velocity_z_north ) / 2
-                                 + velocity_y_north * std::sqrt( this->gamma * pressure_north / u[ north ] )
-                                 * ( std::sqrt( this->gamma * pressure_north / u[ north ] ) - velocity_y_north )
-                                 + this->gamma * pressure_north / u[ north ]
-                                 * ( velocity_y_north - std::sqrt( this->gamma * pressure_north / u[ north ] ) )
-                                 / ( this->gamma - 1 )
-                                 )
+                -hyInverse * ( 
+                                   this->positiveEnergyFlux( density_center, velocity_y_center, velocity_x_center, velocity_z_center, pressure_center)
+                                 - this->positiveEnergyFlux( density_south , velocity_y_south , velocity_x_south , velocity_z_south , pressure_south )
+                                 - this->negativeEnergyFlux( density_center, velocity_y_center, velocity_x_center, velocity_z_center, pressure_center)
+                                 + this->negativeEnergyFlux( density_north , velocity_y_north , velocity_x_north , velocity_z_north , pressure_north ) 
                              ) 
-                -hzInverse * (
-                               ( u[ center ] / ( 2 * this->gamma ) )
-                               * ( 
-                                   ( ( 2 * this->gamma - 1 ) * velocity_z_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                 * ( velocity_x_center * velocity_x_center + velocity_y_center * velocity_y_center + velocity_z_center * velocity_z_center ) / 2
-                                 + velocity_z_center * std::sqrt( this->gamma * pressure_center / u[ center ] )
-                                 * ( velocity_z_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                 + this->gamma * pressure_center / u[ center ] 
-                                 * ( velocity_z_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                 / ( this->gamma - 1 )
-                                 )
-                             - ( u[ down ] / ( 2 * this->gamma ) )
-                               * ( 
-                                   ( ( 2 * this->gamma - 1 ) * velocity_z_down + std::sqrt( this->gamma * pressure_down / u[ down ] ) )
-                                 * ( velocity_x_down * velocity_x_down + velocity_y_down * velocity_y_down + velocity_z_down * velocity_z_down ) / 2
-                                 + velocity_z_down * std::sqrt( this->gamma * pressure_down / u[ down ] )
-                                 * ( velocity_z_down + std::sqrt( this->gamma * pressure_down / u[ down ] ) )
-                                 + this->gamma * pressure_down / u[ down ]
-                                 * ( velocity_z_down + std::sqrt( this->gamma * pressure_down / u[ down ] ) )
-                                 / ( this->gamma - 1 )
-                                 )
-                             - ( u[ center ] / ( 2 * this->gamma ) )
-                               * ( 
-                                   ( velocity_z_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                 * ( velocity_x_center * velocity_x_center + velocity_y_center * velocity_y_center + velocity_z_center * velocity_z_center ) / 2
-                                 + velocity_z_center * std::sqrt( this->gamma * pressure_center / u[ center ] )
-                                 * ( std::sqrt( this->gamma * pressure_center / u[ center ] ) - velocity_z_center )
-                                 + this->gamma * pressure_center / u[ center ]
-                                 * ( velocity_z_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                 / ( this->gamma - 1 )
-                                 )
-                             + ( u[ up ] / ( 2 * this->gamma ) )
-                               * ( 
-                                   ( velocity_z_up - std::sqrt( this->gamma * pressure_up / u[ up ] ) )
-                                 * ( velocity_x_up * velocity_x_up + velocity_y_up * velocity_y_up + velocity_z_up * velocity_z_up ) / 2
-                                 + velocity_z_up * std::sqrt( this->gamma * pressure_up / u[ up ] )
-                                 * ( std::sqrt( this->gamma * pressure_up / u[ up ] ) - velocity_z_up )
-                                 + this->gamma * pressure_up / u[ up ]
-                                 * ( velocity_z_up - std::sqrt( this->gamma * pressure_up / u[ up ] ) )
-                                 / ( this->gamma - 1 )
-                                 )
-                             );
-; 
+                -hyInverse * ( 
+                                   this->positiveEnergyFlux( density_center, velocity_y_center, velocity_x_center, velocity_z_center, pressure_center)
+                                 - this->positiveEnergyFlux( density_down  , velocity_y_down  , velocity_x_down  , velocity_z_down  , pressure_down  )
+                                 - this->negativeEnergyFlux( density_center, velocity_y_center, velocity_x_center, velocity_z_center, pressure_center)
+                                 + this->negativeEnergyFlux( density_up    , velocity_y_up    , velocity_x_up    , velocity_z_up    , pressure_up    ) 
+                             ); 
       }
 
       /*template< typename MeshEntity >
diff --git a/examples/inviscid-flow-sw/UpwindMomentumBase.h b/examples/inviscid-flow-sw/UpwindMomentumBase.h
index 2cb3869d074b8e3eafd53bbf186bb647e6882b6a..5fb86e129fa8ed86217b0f5dfc7c68de800473b1 100644
--- a/examples/inviscid-flow-sw/UpwindMomentumBase.h
+++ b/examples/inviscid-flow-sw/UpwindMomentumBase.h
@@ -46,6 +46,11 @@ class UpwindMomentumBase
       {
           this->velocity = velocity;
       };
+
+      void setDensity( const MeshFunctionPointer& density )
+      {
+          this->density = density;
+      };
       
       void setPressure( const MeshFunctionPointer& pressure )
       {
@@ -85,7 +90,7 @@ class UpwindMomentumBase
          const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
          const RealType& machNumber = velocity_main / speedOfSound;
          if ( machNumber <= -1.0 )
-            return 0;
+            return 0.0;
         else if ( machNumber <= 0.0 )
             return density * speedOfSound / ( 2 * this->gamma ) * ( machNumber + 1.0 ) * velocity_other;
         else if ( machNumber <= 1.0 )
@@ -101,11 +106,11 @@ class UpwindMomentumBase
          if ( machNumber <= -1.0 )
             return density * velocity_main * velocity_other;
         else if ( machNumber <= 0.0 )
-            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( ( 2.0 * this->gamma - 1.0 ) * machNumber - 1.0 ) * velocity_other / speedOfSound;
+            return density * speedOfSound / ( 2 * this->gamma ) * ( ( 2.0 * this->gamma - 1.0 ) * machNumber - 1.0 ) * velocity_other;
         else if ( machNumber <= 1.0 )
-            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( machNumber - 1.0 ) * velocity_other / speedOfSound;
+            return density * speedOfSound / ( 2 * this->gamma ) * ( machNumber - 1.0 ) * velocity_other;
         else 
-            return 0;
+            return 0.0;
       };
 
       protected:
@@ -118,6 +123,8 @@ class UpwindMomentumBase
          
          MeshFunctionPointer pressure;
 
+         MeshFunctionPointer density;
+
 };
 
 } //namespace TNL
diff --git a/examples/inviscid-flow-sw/UpwindMomentumX.h b/examples/inviscid-flow-sw/UpwindMomentumX.h
index 3a6b1abecf3adf5378be5a44d87d4899b6d37d73..ed49dda94585e64f85d820569a757d849757e6ca 100644
--- a/examples/inviscid-flow-sw/UpwindMomentumX.h
+++ b/examples/inviscid-flow-sw/UpwindMomentumX.h
@@ -70,18 +70,24 @@ class UpwindMomentumX< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Ind
          const IndexType& center = entity.getIndex(); 
          const IndexType& east = neighborEntities.template getEntityIndex< 1 >(); 
          const IndexType& west = neighborEntities.template getEntityIndex< -1 >();
+
          const RealType& pressure_center = this->pressure.template getData< DeviceType >()[ center ];
          const RealType& pressure_west = this->pressure.template getData< DeviceType >()[ west ];
          const RealType& pressure_east = this->pressure.template getData< DeviceType >()[ east ];
+
+         const RealType& density_center = this->density.template getData< DeviceType >()[ center ];
+         const RealType& density_west = this->density.template getData< DeviceType >()[ west ];
+         const RealType& density_east = this->density.template getData< DeviceType >()[ east ];
+
          const RealType& velocity_x_center = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ center ];
          const RealType& velocity_x_east = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ];
          const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
          
          return -hxInverse * ( 
-                                 this->positiveMainMomentumFlux( u[ center ], velocity_x_center, pressure_center )
-                               - this->positiveMainMomentumFlux( u[ west   ], velocity_x_west  , pressure_west   )
-                               - this->negativeMainMomentumFlux( u[ center ], velocity_x_center, pressure_center )
-                               + this->negativeMainMomentumFlux( u[ east   ], velocity_x_east  , pressure_east   )
+                                 this->positiveMainMomentumFlux( density_center, velocity_x_center, pressure_center )
+                               - this->positiveMainMomentumFlux( density_west,   velocity_x_west  , pressure_west   )
+                               - this->negativeMainMomentumFlux( density_center, velocity_x_center, pressure_center )
+                               + this->negativeMainMomentumFlux( density_east,   velocity_x_east  , pressure_east   )
                              );
       }
 
@@ -158,6 +164,12 @@ class UpwindMomentumX< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Ind
          const RealType& pressure_east   = this->pressure.template getData< DeviceType >()[ east ];
          const RealType& pressure_north  = this->pressure.template getData< DeviceType >()[ north ];
          const RealType& pressure_south  = this->pressure.template getData< DeviceType >()[ south ];
+         
+         const RealType& density_center = this->density.template getData< DeviceType >()[ center ];
+         const RealType& density_west   = this->density.template getData< DeviceType >()[ west ];
+         const RealType& density_east   = this->density.template getData< DeviceType >()[ east ];
+         const RealType& density_north  = this->density.template getData< DeviceType >()[ north ];
+         const RealType& density_south  = this->density.template getData< DeviceType >()[ south ];
 
          const RealType& velocity_x_center = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ center ];
          const RealType& velocity_x_east   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ];
@@ -166,22 +178,20 @@ class UpwindMomentumX< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Ind
          const RealType& velocity_x_north  = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ north ];
 
          const RealType& velocity_y_center = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ center ];
-         const RealType& velocity_y_east   = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ east ];
-         const RealType& velocity_y_west   = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ west ];
          const RealType& velocity_y_north  = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ north ];
          const RealType& velocity_y_south  = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ];          
          
          return -hxInverse * ( 
-                                 this->positiveMainMomentumFlux( u[ center ], velocity_x_center, pressure_center )
-                               - this->positiveMainMomentumFlux( u[ west   ], velocity_x_west  , pressure_west   )
-                               - this->negativeMainMomentumFlux( u[ center ], velocity_x_center, pressure_center )
-                               + this->negativeMainMomentumFlux( u[ east   ], velocity_x_east  , pressure_east   )
+                                 this->positiveMainMomentumFlux( density_center, velocity_x_center, pressure_center )
+                               - this->positiveMainMomentumFlux( density_west  , velocity_x_west  , pressure_west   )
+                               - this->negativeMainMomentumFlux( density_center, velocity_x_center, pressure_center )
+                               + this->negativeMainMomentumFlux( density_east  , velocity_x_east  , pressure_east   )
                              )
                 -hyInverse * ( 
-                                 this->positiveOtherMomentumFlux( u[ center ], velocity_x_center, velocity_y_center, pressure_center )
-                               - this->positiveOtherMomentumFlux( u[ south  ], velocity_x_south , velocity_y_south , pressure_south  )
-                               - this->negativeOtherMomentumFlux( u[ center ], velocity_x_center, velocity_y_center, pressure_center )
-                               + this->negativeOtherMomentumFlux( u[ north  ], velocity_x_north , velocity_y_north , pressure_north  )
+                                 this->positiveOtherMomentumFlux( density_center, velocity_x_center, velocity_y_center, pressure_center )
+                               - this->positiveOtherMomentumFlux( density_south , velocity_x_south , velocity_y_south , pressure_south  )
+                               - this->negativeOtherMomentumFlux( density_center, velocity_x_center, velocity_y_center, pressure_center )
+                               + this->negativeOtherMomentumFlux( density_north , velocity_x_north , velocity_y_north , pressure_north  )
                              );
 
                                       
@@ -254,9 +264,9 @@ class UpwindMomentumX< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real, Inde
          const IndexType& west   = neighborEntities.template getEntityIndex< -1,  0,  0 >(); 
          const IndexType& north  = neighborEntities.template getEntityIndex<  0,  1,  0 >(); 
          const IndexType& south  = neighborEntities.template getEntityIndex<  0, -1,  0 >();
-         const IndexType& up     = neighborEntities.template getEntityIndex<  0,  0,  1 >(); 
-         const IndexType& down   = neighborEntities.template getEntityIndex<  0,  0, -1 >();
-         
+         const IndexType& up     = neighborEntities.template getEntityIndex<  0,  0,  1 >();
+         const IndexType& down   = neighborEntities.template getEntityIndex<  0,  0, -1 >(); 
+              
          const RealType& pressure_center = this->pressure.template getData< DeviceType >()[ center ];
          const RealType& pressure_west   = this->pressure.template getData< DeviceType >()[ west ];
          const RealType& pressure_east   = this->pressure.template getData< DeviceType >()[ east ];
@@ -265,6 +275,14 @@ class UpwindMomentumX< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real, Inde
          const RealType& pressure_up     = this->pressure.template getData< DeviceType >()[ up ];
          const RealType& pressure_down   = this->pressure.template getData< DeviceType >()[ down ];
          
+         const RealType& density_center = this->density.template getData< DeviceType >()[ center ];
+         const RealType& density_west   = this->density.template getData< DeviceType >()[ west ];
+         const RealType& density_east   = this->density.template getData< DeviceType >()[ east ];
+         const RealType& density_north  = this->density.template getData< DeviceType >()[ north ];
+         const RealType& density_south  = this->density.template getData< DeviceType >()[ south ];
+         const RealType& density_up     = this->density.template getData< DeviceType >()[ up ];
+         const RealType& density_down   = this->density.template getData< DeviceType >()[ down ];
+         
          const RealType& velocity_x_center = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ center ];
          const RealType& velocity_x_east   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ];
          const RealType& velocity_x_west   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
@@ -274,66 +292,30 @@ class UpwindMomentumX< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real, Inde
          const RealType& velocity_x_down   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ down ];
 
          const RealType& velocity_y_center = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ center ];
-         const RealType& velocity_y_east   = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ east ];
-         const RealType& velocity_y_west   = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ west ];
          const RealType& velocity_y_north  = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ north ];
          const RealType& velocity_y_south  = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ];
-         const RealType& velocity_y_up     = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ up ];
-         const RealType& velocity_y_down   = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ down ];
 
          const RealType& velocity_z_center = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ center ];
-         const RealType& velocity_z_east   = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ east ];
-         const RealType& velocity_z_west   = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ west ];
-         const RealType& velocity_z_north  = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ north ];
-         const RealType& velocity_z_south  = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ south ];
          const RealType& velocity_z_up     = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ up ];
          const RealType& velocity_z_down   = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ down ]; 
 
          return -hxInverse * ( 
-                                 u[ center ] / ( 2 * this->gamma )
-                                 * ( 2 * ( this->gamma - 1 ) * velocity_x_center * velocity_x_center
-                                   + ( velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) ) 
-                                   * ( velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                   )
-                               - u[ west   ] / ( 2 * this->gamma )
-                                   * ( 2 * ( this->gamma - 1 ) * velocity_x_west    * velocity_x_west
-                                     + ( velocity_x_west    + std::sqrt( this->gamma * pressure_west   / u[ west ]   ) )
-                                     * ( velocity_x_west    + std::sqrt( this->gamma * pressure_west   / u[ west ]   ) )
-                                     )  
-                                - u[ center ] / ( 2 * this->gamma )
-                                  * ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                  * ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                + u[ east   ] / ( 2 * this->gamma )
-                                  * ( velocity_x_east   - std::sqrt( this->gamma * pressure_center / u[ east   ] ) )
-                                  * ( velocity_x_east   - std::sqrt( this->gamma * pressure_center / u[ east   ] ) )
+                                 this->positiveMainMomentumFlux( density_center, velocity_x_center, pressure_center )
+                               - this->positiveMainMomentumFlux( density_west  , velocity_x_west  , pressure_west   )
+                               - this->negativeMainMomentumFlux( density_center, velocity_x_center, pressure_center )
+                               + this->negativeMainMomentumFlux( density_east  , velocity_x_east  , pressure_east   )
                              )
                 -hyInverse * ( 
-                               u[ center ] / ( 2 * this->gamma )
-                               * ( 2 * ( this->gamma - 1 ) * velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                               * velocity_y_center
-                             - u[ south ] / ( 2 * this->gamma )
-                               * ( 2 * ( this->gamma - 1 ) * velocity_x_south + std::sqrt( this->gamma * pressure_south / u[ south ] ) )
-                               * velocity_y_south
-                             - u[ center ] / ( 2 * this->gamma )
-                               * ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                               * velocity_y_center
-                             + u[ north ] / ( 2 * this->gamma )
-                               * ( velocity_x_north - std::sqrt( this->gamma * pressure_north / u[ north ] ) )
-                               * velocity_y_north
+                                 this->positiveOtherMomentumFlux( density_center, velocity_x_center, velocity_y_center, pressure_center )
+                               - this->positiveOtherMomentumFlux( density_south , velocity_x_south , velocity_y_south , pressure_south  )
+                               - this->negativeOtherMomentumFlux( density_center, velocity_x_center, velocity_y_center, pressure_center )
+                               + this->negativeOtherMomentumFlux( density_north , velocity_x_north , velocity_y_north , pressure_north  )
                              )
                 -hzInverse * ( 
-                               u[ center ] / ( 2 * this->gamma )
-                               * ( 2 * ( this->gamma - 1 ) * velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                               * velocity_z_center
-                             - u[ south ] / ( 2 * this->gamma )
-                               * ( 2 * ( this->gamma - 1 ) * velocity_x_down + std::sqrt( this->gamma * pressure_down / u[ down ] ) )
-                               * velocity_z_south
-                             - u[ center ] / ( 2 * this->gamma )
-                               * ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                               * velocity_z_center
-                             + u[ up ] / ( 2 * this->gamma )
-                               * ( velocity_x_up - std::sqrt( this->gamma * pressure_up / u[ up ] ) )
-                               * velocity_z_up
+                                 this->positiveOtherMomentumFlux( density_center, velocity_x_center, velocity_z_center, pressure_center )
+                               - this->positiveOtherMomentumFlux( density_down  , velocity_x_down  , velocity_z_down  , pressure_down   )
+                               - this->negativeOtherMomentumFlux( density_center, velocity_x_center, velocity_z_center, pressure_center )
+                               + this->negativeOtherMomentumFlux( density_up    , velocity_x_up    , velocity_z_up    , pressure_up     )
                              );
       }
 
diff --git a/examples/inviscid-flow-sw/UpwindMomentumX.h~ b/examples/inviscid-flow-sw/UpwindMomentumX.h~
deleted file mode 100644
index 095ac28c2160b3c9f26b97c8114ccfefbce8b316..0000000000000000000000000000000000000000
--- a/examples/inviscid-flow-sw/UpwindMomentumX.h~
+++ /dev/null
@@ -1,360 +0,0 @@
-/***************************************************************************
-                          UpwindMomentumX.h  -  description
-                             -------------------
-    begin                : Feb 18, 2017
-    copyright            : (C) 2017 by Tomas Oberhuber
-    email                : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/* See Copyright Notice in tnl/Copyright */
-
-
-#pragma once
-
-#include <TNL/Containers/Vector.h>
-#include <TNL/Meshes/Grid.h>
-#include "UpwindMomentumBase.h"
-
-namespace TNL {
-
-template< typename Mesh,
-          typename Real = typename Mesh::RealType,
-          typename Index = typename Mesh::IndexType >
-class UpwindMomentumX
-{
-};
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          typename Index >
-class UpwindMomentumX< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >
-   : public UpwindMomentumBase< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >
-{
-   public:
-
-      typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType;
-      typedef UpwindMomentumBase< MeshType, Real, Index > BaseType;
-      
-      using typename BaseType::RealType;
-      using typename BaseType::IndexType;
-      using typename BaseType::DeviceType;
-      using typename BaseType::CoordinatesType;
-      using typename BaseType::MeshFunctionType;
-      using typename BaseType::MeshFunctionPointer;
-      using typename BaseType::VelocityFieldType;
-      using typename BaseType::VelocityFieldPointer;
-      using BaseType::Dimensions;
-      
-      static String getType()
-      {
-         return String( "UpwindMomentumX< " ) +
-             MeshType::getType() + ", " +
-             TNL::getType< Real >() + ", " +
-             TNL::getType< Index >() + " >"; 
-      }
-      
-
-      template< typename MeshFunction, typename MeshEntity >
-      __cuda_callable__
-      Real operator()( const MeshFunction& u,
-                       const MeshEntity& entity,
-                       const RealType& time = 0.0 ) const
-      {
-         static_assert( MeshEntity::getEntityDimension() == 1, "Wrong mesh entity dimensions." ); 
-         static_assert( MeshFunction::getEntitiesDimension() == 1, "Wrong preimage function" ); 
-         const typename MeshEntity::template NeighborEntities< 1 >& neighborEntities = entity.getNeighborEntities(); 
-
-         const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1 >(); 
-         const IndexType& center = entity.getIndex(); 
-         const IndexType& east = neighborEntities.template getEntityIndex< 1 >(); 
-         const IndexType& west = neighborEntities.template getEntityIndex< -1 >();
-         const RealType& pressure_center = this->pressure.template getData< DeviceType >()[ center ];
-         const RealType& pressure_west = this->pressure.template getData< DeviceType >()[ west ];
-         const RealType& pressure_east = this->pressure.template getData< DeviceType >()[ east ];
-         const RealType& velocity_x_center = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ center ];
-         const RealType& velocity_x_east = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ];
-         const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
-         
-         return -hxInverse * ( 
-                                 this->positiveMainMomentumFlux( u[ center ], velocity_x_center, pressure_center )
-                               - this->positiveMainMomentumFlux( u[ west   ], velocity_x_west  , pressure_west   )
-                               - this->negativeMainMomentumFlux( u[ center ], velocity_x_center, pressure_center )
-                               + this->negativeMainMomentumFlux( u[ east   ], velocity_x_east  , pressure_east   )
-                             );
-      }
-
-      /*template< typename MeshEntity >
-      __cuda_callable__
-      Index getLinearSystemRowLength( const MeshType& mesh,
-                                      const IndexType& index,
-                                      const MeshEntity& entity ) const;
-
-      template< typename MeshEntity, typename Vector, typename MatrixRow >
-      __cuda_callable__
-      void updateLinearSystem( const RealType& time,
-                               const RealType& tau,
-                               const MeshType& mesh,
-                               const IndexType& index,
-                               const MeshEntity& entity,
-                               const MeshFunctionType& u,
-                               Vector& b,
-                               MatrixRow& matrixRow ) const;*/
-};
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          typename Index >
-class UpwindMomentumX< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >
-   : public UpwindMomentumBase< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >
-{
-   public:
-      typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType;
-      typedef UpwindMomentumBase< MeshType, Real, Index > BaseType;
-      
-      using typename BaseType::RealType;
-      using typename BaseType::IndexType;
-      using typename BaseType::DeviceType;
-      using typename BaseType::CoordinatesType;
-      using typename BaseType::MeshFunctionType;
-      using typename BaseType::MeshFunctionPointer;
-      using typename BaseType::VelocityFieldType;
-      using typename BaseType::VelocityFieldPointer;
-      using BaseType::Dimensions;
-      
-      static String getType()
-      {
-         return String( "UpwindMomentumX< " ) +
-             MeshType::getType() + ", " +
-             TNL::getType< Real >() + ", " +
-             TNL::getType< Index >() + " >"; 
-      }      
-
-      template< typename MeshFunction, typename MeshEntity >
-      __cuda_callable__
-      Real operator()( const MeshFunction& u,
-                       const MeshEntity& entity,
-                       const RealType& time = 0.0 ) const
-      {
-         static_assert( MeshEntity::getEntityDimension() == 2, "Wrong mesh entity dimensions." ); 
-         static_assert( MeshFunction::getEntitiesDimension() == 2, "Wrong preimage function" ); 
-         const typename MeshEntity::template NeighborEntities< 2 >& neighborEntities = entity.getNeighborEntities(); 
-
-
-         const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1, 0 >(); 
-         const RealType& hyInverse = entity.getMesh().template getSpaceStepsProducts< 0, -1 >(); 
-
-         const IndexType& center = entity.getIndex(); 
-         const IndexType& east   = neighborEntities.template getEntityIndex<  1,  0 >(); 
-         const IndexType& west   = neighborEntities.template getEntityIndex< -1,  0 >(); 
-         const IndexType& north  = neighborEntities.template getEntityIndex<  0,  1 >(); 
-         const IndexType& south  = neighborEntities.template getEntityIndex<  0, -1 >();
-         
-         const RealType& pressure_center = this->pressure.template getData< DeviceType >()[ center ];
-         const RealType& pressure_west   = this->pressure.template getData< DeviceType >()[ west ];
-         const RealType& pressure_east   = this->pressure.template getData< DeviceType >()[ east ];
-         const RealType& pressure_north  = this->pressure.template getData< DeviceType >()[ north ];
-         const RealType& pressure_south  = this->pressure.template getData< DeviceType >()[ south ];
-
-         const RealType& velocity_x_center = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ center ];
-         const RealType& velocity_x_east   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ];
-         const RealType& velocity_x_west   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
-         const RealType& velocity_x_south  = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ south ];
-         const RealType& velocity_x_north  = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ north ];
-
-         const RealType& velocity_y_center = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ center ];
-         const RealType& velocity_y_east   = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ east ];
-         const RealType& velocity_y_west   = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ west ];
-         const RealType& velocity_y_north  = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ north ];
-         const RealType& velocity_y_south  = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ];          
-         
-         return -hxInverse * ( 
-                                 this->positiveMainMomentumFlux( u[ center ], velocity_x_center, pressure_center )
-                               - this->positiveMainMomentumFlux( u[ west   ], velocity_x_west  , pressure_west   )
-                               - this->negativeMainMomentumFlux( u[ center ], velocity_x_center, pressure_center )
-                               + this->negativeMainMomentumFlux( u[ east   ], velocity_x_east  , pressure_east   )
-                             )
-/*                -hyInverse * ( 
-                                 this->positiveOtherMomentumFlux( u[ center ], velocity_x_center, velocity_y_center, pressure_center )
-                               - this->positiveOtherMomentumFlux( u[ south  ], velocity_x_south , velocity_y_south , pressure_south  )
-                               - this->negativeOtherMomentumFlux( u[ center ], velocity_x_center, velocity_y_center, pressure_center )
-                               + this->negativeOtherMomentumFlux( u[ north  ], velocity_x_north , velocity_y_north , pressure_north  )
-                             )*/;
-
-                                      
-      }
-
-      /*template< typename MeshEntity >
-      __cuda_callable__
-      Index getLinearSystemRowLength( const MeshType& mesh,
-                                      const IndexType& index,
-                                      const MeshEntity& entity ) const;
-
-      template< typename MeshEntity, typename Vector, typename MatrixRow >
-      __cuda_callable__
-      void updateLinearSystem( const RealType& time,
-                               const RealType& tau,
-                               const MeshType& mesh,
-                               const IndexType& index,
-                               const MeshEntity& entity,
-                               const MeshFunctionType& u,
-                               Vector& b,
-                               MatrixRow& matrixRow ) const;*/
-};
-
-template< typename MeshReal,
-          typename Device,
-          typename MeshIndex,
-          typename Real,
-          typename Index >
-class UpwindMomentumX< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real, Index >
-   : public UpwindMomentumBase< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >
-{
-   public:
-      typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;
-      typedef UpwindMomentumBase< MeshType, Real, Index > BaseType;
-      
-      using typename BaseType::RealType;
-      using typename BaseType::IndexType;
-      using typename BaseType::DeviceType;
-      using typename BaseType::CoordinatesType;
-      using typename BaseType::MeshFunctionType;
-      using typename BaseType::MeshFunctionPointer;
-      using typename BaseType::VelocityFieldType;
-      using typename BaseType::VelocityFieldPointer;
-      using BaseType::Dimensions;      
-      
-      static String getType()
-      {
-         return String( "UpwindMomentumX< " ) +
-             MeshType::getType() + ", " +
-             TNL::getType< Real >() + ", " +
-             TNL::getType< Index >() + " >"; 
-      }      
-
-      template< typename MeshFunction, typename MeshEntity >
-      __cuda_callable__
-      Real operator()( const MeshFunction& u,
-                       const MeshEntity& entity,
-                       const RealType& time = 0.0 ) const
-      {
-         static_assert( MeshEntity::getEntityDimension() == 3, "Wrong mesh entity dimensions." ); 
-         static_assert( MeshFunction::getEntitiesDimension() == 3, "Wrong preimage function" ); 
-         const typename MeshEntity::template NeighborEntities< 3 >& neighborEntities = entity.getNeighborEntities(); 
- 
-         const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1, 0,  0 >(); 
-         const RealType& hyInverse = entity.getMesh().template getSpaceStepsProducts< 0, -1,  0 >(); 
-         const RealType& hzInverse = entity.getMesh().template getSpaceStepsProducts< 0,  0, -1 >();
- 
-         const IndexType& center = entity.getIndex(); 
-         const IndexType& east   = neighborEntities.template getEntityIndex<  1,  0,  0 >(); 
-         const IndexType& west   = neighborEntities.template getEntityIndex< -1,  0,  0 >(); 
-         const IndexType& north  = neighborEntities.template getEntityIndex<  0,  1,  0 >(); 
-         const IndexType& south  = neighborEntities.template getEntityIndex<  0, -1,  0 >();
-         const IndexType& up     = neighborEntities.template getEntityIndex<  0,  0,  1 >(); 
-         const IndexType& down   = neighborEntities.template getEntityIndex<  0,  0, -1 >();
-         
-         const RealType& pressure_center = this->pressure.template getData< DeviceType >()[ center ];
-         const RealType& pressure_west   = this->pressure.template getData< DeviceType >()[ west ];
-         const RealType& pressure_east   = this->pressure.template getData< DeviceType >()[ east ];
-         const RealType& pressure_north  = this->pressure.template getData< DeviceType >()[ north ];
-         const RealType& pressure_south  = this->pressure.template getData< DeviceType >()[ south ];
-         const RealType& pressure_up     = this->pressure.template getData< DeviceType >()[ up ];
-         const RealType& pressure_down   = this->pressure.template getData< DeviceType >()[ down ];
-         
-         const RealType& velocity_x_center = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ center ];
-         const RealType& velocity_x_east   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ];
-         const RealType& velocity_x_west   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
-         const RealType& velocity_x_north  = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ north ];
-         const RealType& velocity_x_south  = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ south ];
-         const RealType& velocity_x_up     = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ up ];
-         const RealType& velocity_x_down   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ down ];
-
-         const RealType& velocity_y_center = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ center ];
-         const RealType& velocity_y_east   = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ east ];
-         const RealType& velocity_y_west   = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ west ];
-         const RealType& velocity_y_north  = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ north ];
-         const RealType& velocity_y_south  = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ];
-         const RealType& velocity_y_up     = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ up ];
-         const RealType& velocity_y_down   = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ down ];
-
-         const RealType& velocity_z_center = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ center ];
-         const RealType& velocity_z_east   = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ east ];
-         const RealType& velocity_z_west   = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ west ];
-         const RealType& velocity_z_north  = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ north ];
-         const RealType& velocity_z_south  = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ south ];
-         const RealType& velocity_z_up     = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ up ];
-         const RealType& velocity_z_down   = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ down ]; 
-
-         return -hxInverse * ( 
-                                 u[ center ] / ( 2 * this->gamma )
-                                 * ( 2 * ( this->gamma - 1 ) * velocity_x_center * velocity_x_center
-                                   + ( velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) ) 
-                                   * ( velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                   )
-                               - u[ west   ] / ( 2 * this->gamma )
-                                   * ( 2 * ( this->gamma - 1 ) * velocity_x_west    * velocity_x_west
-                                     + ( velocity_x_west    + std::sqrt( this->gamma * pressure_west   / u[ west ]   ) )
-                                     * ( velocity_x_west    + std::sqrt( this->gamma * pressure_west   / u[ west ]   ) )
-                                     )  
-                                - u[ center ] / ( 2 * this->gamma )
-                                  * ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                  * ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                + u[ east   ] / ( 2 * this->gamma )
-                                  * ( velocity_x_east   - std::sqrt( this->gamma * pressure_center / u[ east   ] ) )
-                                  * ( velocity_x_east   - std::sqrt( this->gamma * pressure_center / u[ east   ] ) )
-                             )
-                -hyInverse * ( 
-                               u[ center ] / ( 2 * this->gamma )
-                               * ( 2 * ( this->gamma - 1 ) * velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                               * velocity_y_center
-                             - u[ south ] / ( 2 * this->gamma )
-                               * ( 2 * ( this->gamma - 1 ) * velocity_x_south + std::sqrt( this->gamma * pressure_south / u[ south ] ) )
-                               * velocity_y_south
-                             - u[ center ] / ( 2 * this->gamma )
-                               * ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                               * velocity_y_center
-                             + u[ north ] / ( 2 * this->gamma )
-                               * ( velocity_x_north - std::sqrt( this->gamma * pressure_north / u[ north ] ) )
-                               * velocity_y_north
-                             )
-                -hzInverse * ( 
-                               u[ center ] / ( 2 * this->gamma )
-                               * ( 2 * ( this->gamma - 1 ) * velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                               * velocity_z_center
-                             - u[ south ] / ( 2 * this->gamma )
-                               * ( 2 * ( this->gamma - 1 ) * velocity_x_down + std::sqrt( this->gamma * pressure_down / u[ down ] ) )
-                               * velocity_z_south
-                             - u[ center ] / ( 2 * this->gamma )
-                               * ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                               * velocity_z_center
-                             + u[ up ] / ( 2 * this->gamma )
-                               * ( velocity_x_up - std::sqrt( this->gamma * pressure_up / u[ up ] ) )
-                               * velocity_z_up
-                             );
-      }
-
-      /*template< typename MeshEntity >
-      __cuda_callable__
-      Index getLinearSystemRowLength( const MeshType& mesh,
-                                      const IndexType& index,
-                                      const MeshEntity& entity ) const;
-
-      template< typename MeshEntity, typename Vector, typename MatrixRow >
-      __cuda_callable__
-      void updateLinearSystem( const RealType& time,
-                               const RealType& tau,
-                               const MeshType& mesh,
-                               const IndexType& index,
-                               const MeshEntity& entity,
-                               const MeshFunctionType& u,
-                               Vector& b,
-                               MatrixRow& matrixRow ) const;*/
-};
-
-
-} // namespace TNL
-
diff --git a/examples/inviscid-flow-sw/UpwindMomentumY.h b/examples/inviscid-flow-sw/UpwindMomentumY.h
index 9859f0618d92d723aba0a4ce3735b00c39e808e5..c2126d43af781289f86999a5f4a7f8d24ad5c6e8 100644
--- a/examples/inviscid-flow-sw/UpwindMomentumY.h
+++ b/examples/inviscid-flow-sw/UpwindMomentumY.h
@@ -142,12 +142,16 @@ class UpwindMomentumY< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Ind
          const RealType& pressure_east   = this->pressure.template getData< DeviceType >()[ east ];
          const RealType& pressure_north  = this->pressure.template getData< DeviceType >()[ north ];
          const RealType& pressure_south  = this->pressure.template getData< DeviceType >()[ south ];
+         
+         const RealType& density_center = this->density.template getData< DeviceType >()[ center ];
+         const RealType& density_west   = this->density.template getData< DeviceType >()[ west ];
+         const RealType& density_east   = this->density.template getData< DeviceType >()[ east ];
+         const RealType& density_north  = this->density.template getData< DeviceType >()[ north ];
+         const RealType& density_south  = this->density.template getData< DeviceType >()[ south ];
 
          const RealType& velocity_x_center = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ center ];
          const RealType& velocity_x_east   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ];
          const RealType& velocity_x_west   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
-         const RealType& velocity_x_south  = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ south ];
-         const RealType& velocity_x_north  = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ north ];
 
          const RealType& velocity_y_center = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ center ];
          const RealType& velocity_y_east   = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ east ];
@@ -156,16 +160,16 @@ class UpwindMomentumY< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Ind
          const RealType& velocity_y_south  = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ];          
          
          return -hxInverse * ( 
-                                 this->positiveOtherMomentumFlux( u[ center ], velocity_y_center, velocity_x_center, pressure_center )
-                               - this->positiveOtherMomentumFlux( u[ west   ], velocity_y_west  , velocity_x_west  , pressure_west   )
-                               - this->negativeOtherMomentumFlux( u[ center ], velocity_y_center, velocity_x_center, pressure_center )
-                               + this->negativeOtherMomentumFlux( u[ east   ], velocity_y_east  , velocity_x_east  , pressure_east   )
+                                 this->positiveOtherMomentumFlux( density_center, velocity_y_center, velocity_x_center, pressure_center )
+                               - this->positiveOtherMomentumFlux( density_west  , velocity_y_west  , velocity_x_west  , pressure_west   )
+                               - this->negativeOtherMomentumFlux( density_center, velocity_y_center, velocity_x_center, pressure_center )
+                               + this->negativeOtherMomentumFlux( density_east  , velocity_y_east  , velocity_x_east  , pressure_east   )
                              )
                 -hyInverse * ( 
-                                 this->positiveMainMomentumFlux( u[ center ], velocity_y_center, pressure_center )
-                               - this->positiveMainMomentumFlux( u[ south  ], velocity_y_south , pressure_south  )
-                               - this->negativeMainMomentumFlux( u[ center ], velocity_y_center, pressure_center )
-                               + this->negativeMainMomentumFlux( u[ north  ], velocity_y_north , pressure_north  )
+                                 this->positiveMainMomentumFlux( density_center, velocity_y_center, pressure_center )
+                               - this->positiveMainMomentumFlux( density_south , velocity_y_south , pressure_south  )
+                               - this->negativeMainMomentumFlux( density_center, velocity_y_center, pressure_center )
+                               + this->negativeMainMomentumFlux( density_north , velocity_y_north , pressure_north  )
                              );
       }
 
@@ -247,13 +251,17 @@ class UpwindMomentumY< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real, Inde
          const RealType& pressure_up     = this->pressure.template getData< DeviceType >()[ up ];
          const RealType& pressure_down   = this->pressure.template getData< DeviceType >()[ down ];
          
+         const RealType& density_center = this->density.template getData< DeviceType >()[ center ];
+         const RealType& density_west   = this->density.template getData< DeviceType >()[ west ];
+         const RealType& density_east   = this->density.template getData< DeviceType >()[ east ];
+         const RealType& density_north  = this->density.template getData< DeviceType >()[ north ];
+         const RealType& density_south  = this->density.template getData< DeviceType >()[ south ];
+         const RealType& density_up     = this->density.template getData< DeviceType >()[ up ];
+         const RealType& density_down   = this->density.template getData< DeviceType >()[ down ];
+         
          const RealType& velocity_x_center = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ center ];
          const RealType& velocity_x_east   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ];
          const RealType& velocity_x_west   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
-         const RealType& velocity_x_north  = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ north ];
-         const RealType& velocity_x_south  = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ south ];
-         const RealType& velocity_x_up     = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ up ];
-         const RealType& velocity_x_down   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ down ];
 
          const RealType& velocity_y_center = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ center ];
          const RealType& velocity_y_east   = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ east ];
@@ -264,58 +272,26 @@ class UpwindMomentumY< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real, Inde
          const RealType& velocity_y_down   = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ down ];
 
          const RealType& velocity_z_center = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ center ];
-         const RealType& velocity_z_east   = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ east ];
-         const RealType& velocity_z_west   = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ west ];
-         const RealType& velocity_z_north  = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ north ];
-         const RealType& velocity_z_south  = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ south ];
          const RealType& velocity_z_up     = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ up ];
          const RealType& velocity_z_down   = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ down ]; 
 
          return -hxInverse * ( 
-                               u[ center ] / ( 2 * this->gamma )
-                               * ( 2 * ( this->gamma - 1 ) * velocity_y_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                               * velocity_x_center
-                             - u[ west ] / ( 2 * this->gamma )
-                               * ( 2 * ( this->gamma - 1 ) * velocity_y_west + std::sqrt( this->gamma * pressure_west / u[ west ] ) )
-                               * velocity_x_west
-                             - u[ center ] / ( 2 * this->gamma )
-                               * ( velocity_y_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                               * velocity_x_center
-                             + u[ east ] / ( 2 * this->gamma )
-                               * ( velocity_y_east - std::sqrt( this->gamma * pressure_east / u[ east ] ) )
-                               * velocity_x_east
+                                 this->positiveOtherMomentumFlux( density_center, velocity_y_center, velocity_x_center, pressure_center )
+                               - this->positiveOtherMomentumFlux( density_west  , velocity_y_west  , velocity_x_west  , pressure_west   )
+                               - this->negativeOtherMomentumFlux( density_center, velocity_y_center, velocity_x_center, pressure_center )
+                               + this->negativeOtherMomentumFlux( density_east  , velocity_y_east  , velocity_x_east  , pressure_east   )
                              )
                 -hyInverse * ( 
-                                 u[ center ] / ( 2 * this->gamma )
-                                 * ( 2 * ( this->gamma - 1 ) * velocity_y_center * velocity_y_center
-                                   + ( velocity_y_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) ) 
-                                   * ( velocity_y_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                   )
-                               - u[ south ] / ( 2 * this->gamma )
-                                   * ( 2 * ( this->gamma - 1 ) * velocity_y_south * velocity_y_south
-                                     + ( velocity_y_south + std::sqrt( this->gamma * pressure_south / u[ south ]   ) )
-                                     * ( velocity_y_south + std::sqrt( this->gamma * pressure_south / u[ south ]   ) )
-                                     )  
-                                - u[ center ] / ( 2 * this->gamma )
-                                  * ( velocity_y_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                  * ( velocity_y_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                + u[ north   ] / ( 2 * this->gamma )
-                                  * ( velocity_y_north   - std::sqrt( this->gamma * pressure_north / u[ north   ] ) )
-                                  * ( velocity_y_north   - std::sqrt( this->gamma * pressure_north / u[ north   ] ) )
+                                 this->positiveMainMomentumFlux( density_center, velocity_y_center, pressure_center )
+                               - this->positiveMainMomentumFlux( density_south , velocity_y_south , pressure_south  )
+                               - this->negativeMainMomentumFlux( density_center, velocity_y_center, pressure_center )
+                               + this->negativeMainMomentumFlux( density_north , velocity_y_north , pressure_north  )
                              )
                 -hzInverse * ( 
-                               u[ center ] / ( 2 * this->gamma )
-                               * ( 2 * ( this->gamma - 1 ) * velocity_y_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                               * velocity_z_center
-                             - u[ south ] / ( 2 * this->gamma )
-                               * ( 2 * ( this->gamma - 1 ) * velocity_y_down + std::sqrt( this->gamma * pressure_down / u[ down ] ) )
-                               * velocity_z_down
-                             - u[ center ] / ( 2 * this->gamma )
-                               * ( velocity_y_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                               * velocity_z_center
-                             + u[ up ] / ( 2 * this->gamma )
-                               * ( velocity_y_up - std::sqrt( this->gamma * pressure_up / u[ up ] ) )
-                               * velocity_z_up
+                                 this->positiveOtherMomentumFlux( density_center, velocity_y_center, velocity_z_center, pressure_center )
+                               - this->positiveOtherMomentumFlux( density_down  , velocity_y_down  , velocity_z_down  , pressure_down   )
+                               - this->negativeOtherMomentumFlux( density_center, velocity_y_center, velocity_z_center, pressure_center )
+                               + this->negativeOtherMomentumFlux( density_up    , velocity_y_up    , velocity_z_up    , pressure_up     )
                              );
       }
 
diff --git a/examples/inviscid-flow-sw/UpwindMomentumZ.h b/examples/inviscid-flow-sw/UpwindMomentumZ.h
index 64842a79f25dadd1be73a2dfb00f166c4885d823..97339e804b3bda5203d0b12feeb59e30249f2327 100644
--- a/examples/inviscid-flow-sw/UpwindMomentumZ.h
+++ b/examples/inviscid-flow-sw/UpwindMomentumZ.h
@@ -208,21 +208,21 @@ class UpwindMomentumZ< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real, Inde
          const RealType& pressure_up     = this->pressure.template getData< DeviceType >()[ up ];
          const RealType& pressure_down   = this->pressure.template getData< DeviceType >()[ down ];
          
+         const RealType& density_center = this->density.template getData< DeviceType >()[ center ];
+         const RealType& density_west   = this->density.template getData< DeviceType >()[ west ];
+         const RealType& density_east   = this->density.template getData< DeviceType >()[ east ];
+         const RealType& density_north  = this->density.template getData< DeviceType >()[ north ];
+         const RealType& density_south  = this->density.template getData< DeviceType >()[ south ];
+         const RealType& density_up     = this->density.template getData< DeviceType >()[ up ];
+         const RealType& density_down   = this->density.template getData< DeviceType >()[ down ];
+         
          const RealType& velocity_x_center = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ center ];
          const RealType& velocity_x_east   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ];
          const RealType& velocity_x_west   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
-         const RealType& velocity_x_north  = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ north ];
-         const RealType& velocity_x_south  = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ south ];
-         const RealType& velocity_x_up     = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ up ];
-         const RealType& velocity_x_down   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ down ];
 
          const RealType& velocity_y_center = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ center ];
-         const RealType& velocity_y_east   = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ east ];
-         const RealType& velocity_y_west   = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ west ];
          const RealType& velocity_y_north  = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ north ];
          const RealType& velocity_y_south  = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ];
-         const RealType& velocity_y_up     = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ up ];
-         const RealType& velocity_y_down   = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ down ];
 
          const RealType& velocity_z_center = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ center ];
          const RealType& velocity_z_east   = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ east ];
@@ -233,50 +233,22 @@ class UpwindMomentumZ< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real, Inde
          const RealType& velocity_z_down   = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ down ]; 
 
          return -hxInverse * ( 
-                               u[ center ] / ( 2 * this->gamma )
-                               * ( ( 2 * this->gamma - 1 ) * velocity_y_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                               * velocity_x_center
-                             - u[ west ] / ( 2 * this->gamma )
-                               * ( ( 2 * this->gamma - 1 ) * velocity_y_west + std::sqrt( this->gamma * pressure_west / u[ west ] ) )
-                               * velocity_x_west
-                             - u[ center ] / ( 2 * this->gamma )
-                               * ( velocity_y_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                               * velocity_x_center
-                             + u[ east ] / ( 2 * this->gamma )
-                               * ( velocity_y_east - std::sqrt( this->gamma * pressure_east / u[ east ] ) )
-                               * velocity_x_east
+                                 this->positiveOtherMomentumFlux( density_center, velocity_z_center, velocity_x_center, pressure_center )
+                               - this->positiveOtherMomentumFlux( density_west  , velocity_z_west  , velocity_x_west  , pressure_west   )
+                               - this->negativeOtherMomentumFlux( density_center, velocity_z_center, velocity_x_center, pressure_center )
+                               + this->negativeOtherMomentumFlux( density_east  , velocity_z_east  , velocity_x_east  , pressure_east   )
                              )
                 -hyInverse * ( 
-                               u[ center ] / ( 2 * this->gamma )
-                               * ( ( 2 * this->gamma - 1 ) * velocity_z_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                               * velocity_y_center
-                             - u[ south ] / ( 2 * this->gamma )
-                               * ( ( 2 * this->gamma - 1 ) * velocity_z_south + std::sqrt( this->gamma * pressure_south / u[ south ] ) )
-                               * velocity_y_south
-                             - u[ center ] / ( 2 * this->gamma )
-                               * ( velocity_z_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                               * velocity_y_center
-                             + u[ north ] / ( 2 * this->gamma )
-                               * ( velocity_z_north - std::sqrt( this->gamma * pressure_north / u[ north ] ) )
-                               * velocity_y_north
+                                 this->positiveOtherMomentumFlux( density_center, velocity_z_center, velocity_y_center, pressure_center )
+                               - this->positiveOtherMomentumFlux( density_south , velocity_z_south , velocity_y_south , pressure_south  )
+                               - this->negativeOtherMomentumFlux( density_center, velocity_z_center, velocity_y_center, pressure_center )
+                               + this->negativeOtherMomentumFlux( density_north , velocity_z_north , velocity_y_north , pressure_north  )
                              )
-                -hyInverse * ( 
-                                 u[ center ] / ( 2 * this->gamma )
-                                 * ( ( 2 * this->gamma - 1 ) * velocity_z_center * velocity_z_center
-                                   + ( velocity_z_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) ) 
-                                   * ( velocity_z_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                   )
-                               - u[ down ] / ( 2 * this->gamma )
-                                   * ( ( 2 * this->gamma - 1 ) * velocity_z_down * velocity_z_down
-                                     + ( velocity_z_down + std::sqrt( this->gamma * pressure_down / u[ down ]   ) )
-                                     * ( velocity_z_down + std::sqrt( this->gamma * pressure_down / u[ down ]   ) )
-                                     )  
-                                - u[ center ] / ( 2 * this->gamma )
-                                  * ( velocity_z_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                  * ( velocity_z_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
-                                + u[ up   ] / ( 2 * this->gamma )
-                                  * ( velocity_z_up   - std::sqrt( this->gamma * pressure_up / u[ up   ] ) )
-                                  * ( velocity_z_up   - std::sqrt( this->gamma * pressure_up / u[ up   ] ) )
+                -hzInverse * ( 
+                                 this->positiveMainMomentumFlux( density_center, velocity_z_center, pressure_center )
+                               - this->positiveMainMomentumFlux( density_down  , velocity_z_down  , pressure_down   )
+                               - this->negativeMainMomentumFlux( density_center, velocity_z_center, pressure_center )
+                               + this->negativeMainMomentumFlux( density_up    , velocity_z_up    , pressure_up     )
                              );
       }
 
diff --git a/examples/inviscid-flow-sw/eulerProblem_impl.h b/examples/inviscid-flow-sw/eulerProblem_impl.h
index e1d5ea9fb1731d293fc43585e19a5e682af6237c..e061360c2cc1ed0a1bcd8783c7a84d8aa64d7828 100644
--- a/examples/inviscid-flow-sw/eulerProblem_impl.h
+++ b/examples/inviscid-flow-sw/eulerProblem_impl.h
@@ -254,6 +254,7 @@ getExplicitUpdate( const RealType& time,
     this->inviscidOperatorsPointer->setTau( tau );
     this->inviscidOperatorsPointer->setVelocity( this->velocity );
     this->inviscidOperatorsPointer->setPressure( this->pressure );
+    this->inviscidOperatorsPointer->setDensity( this->conservativeVariables->getDensity() );
     this->inviscidOperatorsPointer->setGamma( this->gamma );
 
 //   this->pressure->write( "pressure2", "gnuplot" );
@@ -277,7 +278,7 @@ getExplicitUpdate( const RealType& time,
    explicitUpdaterMomentumX.setBoundaryConditions( this->boundaryConditionPointer );
    explicitUpdaterMomentumX.setRightHandSide( this->rightHandSidePointer );   
    explicitUpdaterMomentumX.template update< typename Mesh::Cell >( time, tau, mesh,
-                                                              this->conservativeVariables->getDensity(), // uRhoVelocityX,
+                                                           ( *this->conservativeVariables->getMomentum() )[ 0 ], // uRhoVelocityX,
                                                            ( *this->conservativeVariablesRHS->getMomentum() )[ 0 ] ); //, fuRhoVelocityX );
 
    if( Dimensions > 1 )
@@ -287,7 +288,7 @@ getExplicitUpdate( const RealType& time,
       explicitUpdaterMomentumY.setBoundaryConditions( this->boundaryConditionPointer );
       explicitUpdaterMomentumY.setRightHandSide( this->rightHandSidePointer );         
       explicitUpdaterMomentumY.template update< typename Mesh::Cell >( time, tau, mesh,
-                                                                 this->conservativeVariables->getDensity(), // uRhoVelocityX,
+                                                              ( *this->conservativeVariables->getMomentum() )[ 1 ], // uRhoVelocityX,
                                                               ( *this->conservativeVariablesRHS->getMomentum() )[ 1 ] ); //, fuRhoVelocityX );
    }
    
@@ -298,7 +299,7 @@ getExplicitUpdate( const RealType& time,
       explicitUpdaterMomentumZ.setBoundaryConditions( this->boundaryConditionPointer );
       explicitUpdaterMomentumZ.setRightHandSide( this->rightHandSidePointer );               
       explicitUpdaterMomentumZ.template update< typename Mesh::Cell >( time, tau, mesh,
-                                                                 this->conservativeVariables->getDensity(), // uRhoVelocityX,
+                                                              ( *this->conservativeVariables->getMomentum() )[ 2 ], // uRhoVelocityX,
                                                               ( *this->conservativeVariablesRHS->getMomentum() )[ 2 ] ); //, fuRhoVelocityX );
    }
    
@@ -311,7 +312,7 @@ getExplicitUpdate( const RealType& time,
    explicitUpdaterEnergy.setBoundaryConditions( this->boundaryConditionPointer );
    explicitUpdaterEnergy.setRightHandSide( this->rightHandSidePointer );                  
    explicitUpdaterEnergy.template update< typename Mesh::Cell >( time, tau, mesh,
-                                                           this->conservativeVariables->getDensity(), // uRhoVelocityX,
+                                                           this->conservativeVariablesRHS->getEnergy(), // uRhoVelocityX,
                                                            this->conservativeVariablesRHS->getEnergy() ); //, fuRhoVelocityX );
 
 /*   this->pressure->write( "pressure3", "gnuplot" );