diff --git a/examples/inviscid-flow-sw/UpwindContinuity.h b/examples/inviscid-flow-sw/UpwindContinuity.h
index 4439d831381b0b7c7b8b09c826f481ba6c376625..a375f8b9597f735873a9a148e6fb3a3aa9244830 100644
--- a/examples/inviscid-flow-sw/UpwindContinuity.h
+++ b/examples/inviscid-flow-sw/UpwindContinuity.h
@@ -65,7 +65,7 @@ class UpwindContinuityBase
           this->velocity = velocity;
       };
 
-      const RealType& positiveDensityFlux( const RealType& density, const RealType& velocity, const RealType& pressure )
+      RealType positiveDensityFlux( const RealType& density, const RealType& velocity, const RealType& pressure ) const
       {
          const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
          const RealType& machNumber = velocity / speedOfSound;
@@ -77,9 +77,9 @@ class UpwindContinuityBase
             return density * speedOfSound / ( 2 * this->gamma ) * ( ( 2.0 * this->gamma - 1.0 ) * machNumber + 1.0 );
         else 
             return density * velocity;
-      }
+      };
 
-      const RealType& negativeDensityFlux( const RealType& density, const RealType& velocity, const RealType& pressure )
+      RealType negativeDensityFlux( const RealType& density, const RealType& velocity, const RealType& pressure ) const
       {
          const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
          const RealType& machNumber = velocity / speedOfSound;
@@ -91,7 +91,12 @@ class UpwindContinuityBase
             return density * speedOfSound / ( 2 * this->gamma ) * ( machNumber - 1.0 );
         else 
             return 0.0;
-      }
+      };
+      
+      RealType multiply (const RealType& a, const RealType& b ) const
+      {
+         return a * b;
+      };
       
 
       protected:
@@ -162,11 +167,11 @@ class UpwindContinuity< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, In
          const RealType& velocity_x_east   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ];
 
          return -hxInverse * (
-                                   UpwindContinuity::positiveDensityFlux( u[ center ], velocity_x_center, pressure_center )
-                                -  UpwindContinuity::positiveDensityFlux( u[ west   ], velocity_x_west  , pressure_west   )
-                                -  UpwindContinuity::negativeDensityFlux( u[ center ], velocity_x_center, pressure_center )
-                                +  UpwindContinuity::negativeDensityFlux( u[ east   ], velocity_x_east  , pressure_east   )
-//                                  u[ center ] / ( 2 * this->gamma ) * ( ( 2 * this->gamma - 1 ) * velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
+                                   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   )
+//                                  u[ center ] / ( 2 * this->gamma ) * ( ( 2 * this->gamma - 1 ) * velocity_x_center + std::sqrt( this->multiply( 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   ] ) )
diff --git a/examples/inviscid-flow-sw/UpwindEnergy.h b/examples/inviscid-flow-sw/UpwindEnergy.h
index 8029f3dc6999c05a3c284800825eb587f6a8c910..7ecc6449c2c004bbc30853269abfc1d1265d200d 100644
--- a/examples/inviscid-flow-sw/UpwindEnergy.h
+++ b/examples/inviscid-flow-sw/UpwindEnergy.h
@@ -67,7 +67,7 @@ class UpwindEnergyBase
       void setArtificialViscosity( const RealType& artificialViscosity )
       {
          this->artificialViscosity = artificialViscosity;
-      }
+      };
 
       const RealType& positiveEnergyFlux( const RealType& density, const RealType& velocity_main, const RealType& velocity_other1, const RealType& velocity_other2, const RealType& pressure )
       {
@@ -78,44 +78,44 @@ class UpwindEnergyBase
         else if ( machNumber <= 0.0 )
             return density * speedOfSound / ( 2 * this->gamma ) 
                  * ( 
-                     ( speedOfSound * ( machNumber + 1.0 ) * speedOfSound * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) );
-                   + ( speedOfSound * speedOfSound * machNumber * ( - 1.0 - machNumber )
+                     ( speedOfSound * ( machNumber + 1.0 ) * speedOfSound * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) )
+                   + ( speedOfSound * speedOfSound * machNumber * ( - 1.0 - machNumber ) )
                    + ( speedOfSound * speedOfSound * ( machNumber + 1.0 ) / ( this->gamma - 1.0 ) )
                    );
         else if ( machNumber <= 1.0 )
             return density * speedOfSound / ( 2 * this->gamma ) 
                  * ( 
-                     ( speedOfSound * ( ( 2.0 * this->gamma - 1.0 ) * machNumber + 1.0 ) ) * speedOfSound * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) );
-                   + ( speedOfSound * speedOfSound * machNumber * ( machNumber + 1.0 )
+                     ( speedOfSound * ( ( 2.0 * this->gamma - 1.0 ) * machNumber + 1.0 ) * speedOfSound * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) )
+                   + ( speedOfSound * speedOfSound * machNumber * ( machNumber + 1.0 ) )
                    + ( speedOfSound * speedOfSound * ( machNumber + 1.0 ) / ( this->gamma - 1.0 ) )
                    );
         else   
-            return ( pressure / ( this->gamma - 1.0 ) + 0.5 * density * ( velocity_main * velocity_main + velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 );
-      }
+            return ( pressure / ( this->gamma - 1.0 ) + 0.5 * density * ( velocity_main * velocity_main + velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) );
+      };
 
-      const RealType& negativeEnergyFlux( const RealType& density, const RealType& velocity, const RealType& pressure )
+      const RealType& negativeEnergyFlux( const RealType& density, const RealType& velocity_main, const RealType& velocity_other1, const RealType& velocity_other2, const RealType& pressure )
       {
          const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
          const RealType& machNumber = velocity / speedOfSound;
          if ( machNumber <= -1.0 )
-            return ( pressure / ( this->gamma - 1.0 ) + 0.5 * density * ( velocity_main * velocity_main + velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 );
+            return ( 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 / ( 2 * this->gamma ) 
                  * ( 
-                     ( speedOfSound * ( ( 2.0 * this->gamma - 1.0 ) * machNumber - 1.0 )  * speedOfSound * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) );
-                   + ( speedOfSound * speedOfSound * machNumber * ( machNumber - 1.0 )
+                     ( speedOfSound * ( ( 2.0 * this->gamma - 1.0 ) * machNumber - 1.0 )  * speedOfSound * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) )
+                   + ( speedOfSound * speedOfSound * machNumber * ( machNumber - 1.0 ) )
                    + ( speedOfSound * speedOfSound * ( machNumber - 1.0 ) / ( this->gamma - 1.0 ) )
-                   )
+                   );
         else if ( machNumber <= 1.0 )
             return density * speedOfSound / ( 2 * this->gamma ) 
                  * ( 
-                     ( speedOfSound * ( machNumber - 1.0 ) * speedOfSound * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) );
-                   + ( speedOfSound * speedOfSound * machNumber * ( 1.0 - machNumber )
+                     ( speedOfSound * ( machNumber - 1.0 ) * speedOfSound * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) )
+                   + ( speedOfSound * speedOfSound * machNumber * ( 1.0 - machNumber ) )
                    + ( speedOfSound * speedOfSound * ( machNumber - 1.0 ) / ( this->gamma - 1.0 ) )
                    );
         else 
             return 0.0;
-      }      
+      };      
 
       protected:
          
diff --git a/examples/inviscid-flow-sw/UpwindMomentumBase.h b/examples/inviscid-flow-sw/UpwindMomentumBase.h
index dd33e3a9783bbb95b9a2530582f63f8730295f37..c9bb5745fce068c1f3bc0aa1b66d367d76b99c5a 100644
--- a/examples/inviscid-flow-sw/UpwindMomentumBase.h
+++ b/examples/inviscid-flow-sw/UpwindMomentumBase.h
@@ -59,26 +59,26 @@ class UpwindMomentumBase
          if ( machNumber <= -1.0 )
             return 0;
         else if ( machNumber <= 0.0 )
-            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( ( machNumber + 1.0 ) * machNumber + ( - machNumber - 1.0 );
+            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( ( machNumber + 1.0 ) * machNumber + ( - machNumber - 1.0 ) );
         else if ( machNumber <= 1.0 )
-            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( ( ( 2.0 * this->gamma - 1.0 ) * machNumber + 1.0 ) * machNumber + (machNumber + 1.0 );
+            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( ( ( 2.0 * this->gamma - 1.0 ) * machNumber + 1.0 ) * machNumber + (machNumber + 1.0 ) );
         else 
             return density * velocity * velocity + pressure;
-      }
+      };
 
       const RealType& negativeMainMomentumFlux( const RealType& density, const RealType& velocity, const RealType& pressure )
       {
          const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
          const RealType& machNumber = velocity / speedOfSound;
          if ( machNumber <= -1.0 )
-            return density * velocity;
+            return density * velocity * velocity + pressure;
         else if ( machNumber <= 0.0 )
-            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( ( ( 2.0 * this->gamma - 1.0 ) * machNumber - 1.0 ) * machNumber + (machNumber - 1.0 );
+            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( ( ( 2.0 * this->gamma - 1.0 ) * machNumber - 1.0 ) * machNumber + (machNumber - 1.0 ) );
         else if ( machNumber <= 1.0 )
-            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( ( machNumber - 1.0 ) * machNumber + ( - machNumber + 1.0 );
+            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( ( machNumber - 1.0 ) * machNumber + ( - machNumber + 1.0 ) );
         else 
-            return density * velocity * velocity + pressure;
-      }
+            return 0; 
+      };
 
       const RealType& positiveOtherMomentumFlux( const RealType& density, const RealType& velocity_main, const RealType& velocity_other, const RealType& pressure )
       {
@@ -92,7 +92,7 @@ class UpwindMomentumBase
             return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( ( 2.0 * this->gamma - 1.0 ) * machNumber + 1.0 ) * velocity_other / speedOfSound;
         else 
             return density * velocity_main * velocity_other;
-      }
+      };
 
       const RealType& negativeOtherMomentumFlux( const RealType& density, const RealType& velocity_main, const RealType& velocity_other, const RealType& pressure )
       {
@@ -101,12 +101,12 @@ 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 * speedOfSound / ( 2 * this->gamma ) * ( ( 2.0 * this->gamma - 1.0 ) * machNumber - 1.0 ) * velocity_other / speedOfSound;
         else if ( machNumber <= 1.0 )
-            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( machNumber + 1.0 ) * velocity_other / speedOfSound;
+            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( machNumber - 1.0 ) * velocity_other / speedOfSound;
         else 
             return 0;
-      }
+      };
 
       protected: