diff --git a/examples/inviscid-flow-sw/RiemannProblemInitialCondition.h b/examples/inviscid-flow-sw/RiemannProblemInitialCondition.h
index bfa458a8bb5cdecf9a1b1013f332207ebc3db66e..85fbc8af01994495e2c1f2f7d95f8155ee216684 100644
--- a/examples/inviscid-flow-sw/RiemannProblemInitialCondition.h
+++ b/examples/inviscid-flow-sw/RiemannProblemInitialCondition.h
@@ -949,7 +949,7 @@ class RiemannProblemInitialCondition
0.0, 0.0, 0.0, //double preNEUVelocityX, double preNEUVelocityY,double preNEUVelocityZ,
0.0, 0.0, 0.0, //double preSEUVelocityX, double preSEUVelocityY,double preSEUVelocityZ,
0.0, 0.0, 0.0, //double preNEDVelocityX, double preNEDVelocityY,double preNEDVelocityZ,
- -19.57945, 0.0, 0.0 //double preSEDVelocityX, double preSEDVelocityY,double preSEDVelocityZ
+ -19.59745, 0.0, 0.0 //double preSEDVelocityX, double preSEDVelocityY,double preSEDVelocityZ
);
if(initial == prefix + "1D_4")
predefinedInitialCondition( 1.666, 0.4, 0.0, 0.0, // double preGamma, double preDiscX, double preDiscY, double preDiscZ,
diff --git a/examples/inviscid-flow-sw/UpwindEnergy.h b/examples/inviscid-flow-sw/UpwindEnergy.h
index 7ecc6449c2c004bbc30853269abfc1d1265d200d..29d74e32e94c46e2e5760195a4c910731985cf18 100644
--- a/examples/inviscid-flow-sw/UpwindEnergy.h
+++ b/examples/inviscid-flow-sw/UpwindEnergy.h
@@ -69,49 +69,49 @@ class UpwindEnergyBase
this->artificialViscosity = artificialViscosity;
};
- const RealType& positiveEnergyFlux( const RealType& density, const RealType& velocity_main, const RealType& velocity_other1, const RealType& velocity_other2, const RealType& pressure )
+ 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 / ( 2 * this->gamma )
+ return density * speedOfSound * speedOfSound * speedOfSound / ( 2.0 * 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 * speedOfSound * ( machNumber + 1.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 / ( 2 * this->gamma )
+ return density * speedOfSound * speedOfSound * speedOfSound / ( 2.0 * 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 * speedOfSound * ( machNumber + 1.0 ) / ( this->gamma - 1.0 ) )
+ ( ( ( 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 ( pressure / ( this->gamma - 1.0 ) + 0.5 * density * ( velocity_main * velocity_main + velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) );
+ return velocity_main * ( pressure + 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_main, const RealType& velocity_other1, const RealType& velocity_other2, const RealType& pressure )
+ 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 / speedOfSound;
+ const RealType& machNumber = velocity_main / 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 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 / ( 2 * this->gamma )
+ return density * speedOfSound * speedOfSound * speedOfSound / ( 2.0 * 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 * speedOfSound * ( machNumber - 1.0 ) / ( this->gamma - 1.0 ) )
+ ( ( ( 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 / ( 2 * this->gamma )
+ return density * speedOfSound * speedOfSound * speedOfSound / ( 2.0 * 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 * speedOfSound * ( machNumber - 1.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 0.0;
@@ -185,52 +185,10 @@ class UpwindEnergy< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index
const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
return -hxInverse * (
- ( u[ center ] / ( 2 * this->gamma ) )
- * (
- ( this->gamma - 1 ) * velocity_x_center * 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 ] ) )
- * ( velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
- / 2
- + ( 3 - this->gamma )
- / ( 2 * ( this->gamma - 1 ) )
- * ( velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
- * ( this->gamma * pressure_center / u[ center ] )
- )
- - ( u[ west ] / ( 2 * this->gamma ) )
- * (
- ( this->gamma - 1 ) * velocity_x_west * 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 ] ) )
- * ( velocity_x_west + std::sqrt( this->gamma * pressure_west / u[ west ] ) )
- / 2
- + ( 3 - this->gamma )
- / ( 2 * ( this->gamma - 1 ) )
- * ( velocity_x_west + std::sqrt( this->gamma * pressure_west / u[ west ] ) )
- * ( 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 ] ) )
- * ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
- / 2
- + ( 3 - this->gamma )
- / ( this->gamma - 1 )
- * ( this->gamma * pressure_center / u[ center ] )
- / 2
- )
- + u[ east ] / ( 2 * this->gamma )
- * ( velocity_x_east - std::sqrt( this->gamma * pressure_east / u[ east ] ) )
- * (
- ( velocity_x_east - std::sqrt( this->gamma * pressure_east / u[ east ] ) )
- * ( velocity_x_east - std::sqrt( this->gamma * pressure_east / u[ east ] ) )
- / 2
- + ( 3 - this->gamma )
- / ( this->gamma - 1 )
- * ( this->gamma * pressure_east / u[ east ] )
- / 2
- )*/
+ 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 )
);
}
diff --git a/examples/inviscid-flow-sw/UpwindMomentumBase.h b/examples/inviscid-flow-sw/UpwindMomentumBase.h
index c9bb5745fce068c1f3bc0aa1b66d367d76b99c5a..29d4612c1df34e60f9bf08771673cbde2930c474 100644
--- a/examples/inviscid-flow-sw/UpwindMomentumBase.h
+++ b/examples/inviscid-flow-sw/UpwindMomentumBase.h
@@ -52,35 +52,35 @@ class UpwindMomentumBase
this->pressure = pressure;
};
- const RealType& positiveMainMomentumFlux( const RealType& density, const RealType& velocity, const RealType& pressure )
+ RealType positiveMainMomentumFlux( const RealType& density, const RealType& velocity, const RealType& pressure ) const
{
const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
const RealType& machNumber = velocity / speedOfSound;
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 + 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 * machNumber + (machNumber + 1.0 ) * (machNumber + 1.0 ) );
else
return density * velocity * velocity + pressure;
};
- const RealType& negativeMainMomentumFlux( const RealType& density, const RealType& velocity, const RealType& pressure )
+ RealType negativeMainMomentumFlux( const RealType& density, const RealType& velocity, const RealType& pressure ) const
{
const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
const RealType& machNumber = velocity / speedOfSound;
if ( machNumber <= -1.0 )
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 * machNumber + (machNumber - 1.0 ) * (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 - 1.0 ) );
else
return 0;
};
- const RealType& positiveOtherMomentumFlux( const RealType& density, const RealType& velocity_main, const RealType& velocity_other, const RealType& pressure )
+ RealType positiveOtherMomentumFlux( const RealType& density, const RealType& velocity_main, const RealType& velocity_other, const RealType& pressure ) const
{
const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
const RealType& machNumber = velocity_main / speedOfSound;
@@ -94,7 +94,7 @@ class UpwindMomentumBase
return density * velocity_main * velocity_other;
};
- const RealType& negativeOtherMomentumFlux( const RealType& density, const RealType& velocity_main, const RealType& velocity_other, const RealType& pressure )
+ RealType negativeOtherMomentumFlux( const RealType& density, const RealType& velocity_main, const RealType& velocity_other, const RealType& pressure ) const
{
const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
const RealType& machNumber = velocity_main / speedOfSound;
diff --git a/examples/inviscid-flow-sw/UpwindMomentumX.h b/examples/inviscid-flow-sw/UpwindMomentumX.h
index 662e0707f4b2846da8206d34c59a462503fe3e5a..c9cafec4413c6484e41d0968d77bdda463b00c5b 100644
--- a/examples/inviscid-flow-sw/UpwindMomentumX.h
+++ b/examples/inviscid-flow-sw/UpwindMomentumX.h
@@ -78,24 +78,10 @@ class UpwindMomentumX< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Ind
const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
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_east / u[ east ] ) )
- * ( velocity_x_east - std::sqrt( this->gamma * pressure_east / u[ east ] ) )*/
+ 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 )
);
}
diff --git a/examples/inviscid-flow-sw/UpwindMomentumX.h~ b/examples/inviscid-flow-sw/UpwindMomentumX.h~
index 4b941f23ae47b40dabd3d88dfd58c0d6d90eb039..7f72b31daefabb61fe570b56fc7c41b02cbb5f22 100644
--- a/examples/inviscid-flow-sw/UpwindMomentumX.h~
+++ b/examples/inviscid-flow-sw/UpwindMomentumX.h~
@@ -78,24 +78,10 @@ class UpwindMomentumX< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Ind
const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
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_east / u[ east ] ) )
- * ( velocity_x_east - std::sqrt( this->gamma * pressure_east / u[ east ] ) )
+ this->positiveMainMomentumFlux( u[ center ], velocity_x_center, pressure_center )
+ - this->positiveMainMomentumFlux( u[ west ], velocity_x_west , pressure_west )
+ - this->positiveMainMomentumFlux( u[ center ], velocity_x_center, pressure_center )
+ + this->positiveMainMomentumFlux( u[ east ], velocity_x_east , pressure_east )
);
}
diff --git a/examples/inviscid-flow/RiemannProblemInitialCondition.h b/examples/inviscid-flow/RiemannProblemInitialCondition.h
index a036747f255fd1477c02a55366d2f57ef2f7f990..b1819b1eb86c014f1689945b1dcac44b6cbdbfcd 100644
--- a/examples/inviscid-flow/RiemannProblemInitialCondition.h
+++ b/examples/inviscid-flow/RiemannProblemInitialCondition.h
@@ -949,7 +949,7 @@ class RiemannProblemInitialCondition
0.0, 0.0, 0.0, //double preNEUVelocityX, double preNEUVelocityY,double preNEUVelocityZ,
0.0, 0.0, 0.0, //double preSEUVelocityX, double preSEUVelocityY,double preSEUVelocityZ,
0.0, 0.0, 0.0, //double preNEDVelocityX, double preNEDVelocityY,double preNEDVelocityZ,
- -19.57945, 0.0, 0.0 //double preSEDVelocityX, double preSEDVelocityY,double preSEDVelocityZ
+ -19.59745, 0.0, 0.0 //double preSEDVelocityX, double preSEDVelocityY,double preSEDVelocityZ
);
if(initial == prefix + "1D_4")
predefinedInitialCondition( 1.666, 0.4, 0.0, 0.0, // double preGamma, double preDiscX, double preDiscY, double preDiscZ,