Skip to content
Snippets Groups Projects
Commit 8286beec authored by Jan Schäfer's avatar Jan Schäfer
Browse files

plne funkcni 1D, odzkouseno

parent bd1ccb18
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 43 additions and 113 deletions
examples/inviscid-flow-sw/RiemannProblemInitialCondition.h
+ 1
− 1
View file @ 8286beec
...@@ -949,7 +949,7 @@ class RiemannProblemInitialCondition ...@@ -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 preNEUVelocityX, double preNEUVelocityY,double preNEUVelocityZ,
0.0, 0.0, 0.0, //double preSEUVelocityX, double preSEUVelocityY,double preSEUVelocityZ, 0.0, 0.0, 0.0, //double preSEUVelocityX, double preSEUVelocityY,double preSEUVelocityZ,
0.0, 0.0, 0.0, //double preNEDVelocityX, double preNEDVelocityY,double preNEDVelocityZ, 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") if(initial == prefix + "1D_4")
predefinedInitialCondition( 1.666, 0.4, 0.0, 0.0, // double preGamma, double preDiscX, double preDiscY, double preDiscZ, predefinedInitialCondition( 1.666, 0.4, 0.0, 0.0, // double preGamma, double preDiscX, double preDiscY, double preDiscZ,
......
examples/inviscid-flow-sw/UpwindEnergy.h
+ 25
− 67
View file @ 8286beec
...@@ -69,49 +69,49 @@ class UpwindEnergyBase ...@@ -69,49 +69,49 @@ class UpwindEnergyBase
this->artificialViscosity = 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 ) 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& speedOfSound = std::sqrt( this->gamma * pressure / density );
const RealType& machNumber = velocity_main / speedOfSound; const RealType& machNumber = velocity_main / speedOfSound;
if ( machNumber <= -1.0 ) if ( machNumber <= -1.0 )
return 0.0; return 0.0;
else if ( machNumber <= 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 ) ) ) ( ( machNumber + 1.0 ) * 0.5 * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) )
+ ( speedOfSound * speedOfSound * machNumber * ( - 1.0 - machNumber ) ) + ( machNumber * ( machNumber + 1.0 ) )
+ ( speedOfSound * speedOfSound * ( machNumber + 1.0 ) / ( this->gamma - 1.0 ) ) + ( ( machNumber + 1.0 ) / ( this->gamma - 1.0 ) )
); );
else if ( machNumber <= 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 ) ) ) ( ( ( 2.0 * this->gamma - 1.0 ) * machNumber + 1.0 ) * 0.5 * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) )
+ ( speedOfSound * speedOfSound * machNumber * ( machNumber + 1.0 ) ) + ( machNumber * ( machNumber + 1.0 ) )
+ ( speedOfSound * speedOfSound * ( machNumber + 1.0 ) / ( this->gamma - 1.0 ) ) + ( ( machNumber + 1.0 ) / ( this->gamma - 1.0 ) )
); );
else 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& speedOfSound = std::sqrt( this->gamma * pressure / density );
const RealType& machNumber = velocity / speedOfSound; const RealType& machNumber = velocity_main / speedOfSound;
if ( machNumber <= -1.0 ) 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 ) 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 ) ) ) ( ( ( 2.0 * this->gamma - 1.0 ) * machNumber - 1.0 ) * 0.5 * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) )
+ ( speedOfSound * speedOfSound * machNumber * ( machNumber - 1.0 ) ) - ( machNumber * ( machNumber - 1.0 ) )
+ ( speedOfSound * speedOfSound * ( machNumber - 1.0 ) / ( this->gamma - 1.0 ) ) + ( ( machNumber - 1.0 ) / ( this->gamma - 1.0 ) )
); );
else if ( machNumber <= 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 ) ) ) ( ( machNumber - 1.0 ) * 0.5 * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) )
+ ( speedOfSound * speedOfSound * machNumber * ( 1.0 - machNumber ) ) - ( machNumber * ( machNumber - 1.0 ) )
+ ( speedOfSound * speedOfSound * ( machNumber - 1.0 ) / ( this->gamma - 1.0 ) ) + ( ( machNumber - 1.0 ) / ( this->gamma - 1.0 ) )
); );
else else
return 0.0; return 0.0;
...@@ -185,52 +185,10 @@ class UpwindEnergy< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index ...@@ -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 ]; const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
return -hxInverse * ( return -hxInverse * (
( u[ center ] / ( 2 * this->gamma ) ) this->positiveEnergyFlux( u[ center ], velocity_x_center, 0, 0, pressure_center)
* ( - this->positiveEnergyFlux( u[ west ], velocity_x_west , 0, 0, pressure_west )
( this->gamma - 1 ) * velocity_x_center * velocity_x_center * velocity_x_center - this->negativeEnergyFlux( u[ center ], velocity_x_center, 0, 0, pressure_center)
+ ( velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) ) + this->negativeEnergyFlux( u[ east ], velocity_x_east , 0, 0, pressure_east )
* ( 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
)*/
); );
} }
......
examples/inviscid-flow-sw/UpwindMomentumBase.h
+ 8
− 8
View file @ 8286beec
...@@ -52,35 +52,35 @@ class UpwindMomentumBase ...@@ -52,35 +52,35 @@ class UpwindMomentumBase
this->pressure = pressure; 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& speedOfSound = std::sqrt( this->gamma * pressure / density );
const RealType& machNumber = velocity / speedOfSound; const RealType& machNumber = velocity / speedOfSound;
if ( machNumber <= -1.0 ) if ( machNumber <= -1.0 )
return 0; return 0;
else if ( machNumber <= 0.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 ) 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 else
return density * velocity * velocity + pressure; 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& speedOfSound = std::sqrt( this->gamma * pressure / density );
const RealType& machNumber = velocity / speedOfSound; const RealType& machNumber = velocity / speedOfSound;
if ( machNumber <= -1.0 ) if ( machNumber <= -1.0 )
return density * velocity * velocity + pressure; return density * velocity * velocity + pressure;
else if ( machNumber <= 0.0 ) 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 ) 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 else
return 0; 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& speedOfSound = std::sqrt( this->gamma * pressure / density );
const RealType& machNumber = velocity_main / speedOfSound; const RealType& machNumber = velocity_main / speedOfSound;
...@@ -94,7 +94,7 @@ class UpwindMomentumBase ...@@ -94,7 +94,7 @@ class UpwindMomentumBase
return density * velocity_main * velocity_other; 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& speedOfSound = std::sqrt( this->gamma * pressure / density );
const RealType& machNumber = velocity_main / speedOfSound; const RealType& machNumber = velocity_main / speedOfSound;
......
examples/inviscid-flow-sw/UpwindMomentumX.h
+ 4
− 18
View file @ 8286beec
...@@ -78,24 +78,10 @@ class UpwindMomentumX< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Ind ...@@ -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 ]; const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
return -hxInverse * ( return -hxInverse * (
( u[ center ] / ( 2 * this->gamma ) ) this->positiveMainMomentumFlux( u[ center ], velocity_x_center, pressure_center )
* ( - this->positiveMainMomentumFlux( u[ west ], velocity_x_west , pressure_west )
2 * ( this->gamma - 1 ) * velocity_x_center * velocity_x_center - this->negativeMainMomentumFlux( u[ center ], velocity_x_center, pressure_center )
+ ( velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) ) + this->negativeMainMomentumFlux( u[ east ], velocity_x_east , pressure_east )
* ( 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 ] ) )*/
); );
} }
......
examples/inviscid-flow-sw/UpwindMomentumX.h~
+ 4
− 18
View file @ 8286beec
...@@ -78,24 +78,10 @@ class UpwindMomentumX< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Ind ...@@ -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 ]; const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
return -hxInverse * ( return -hxInverse * (
( u[ center ] / ( 2 * this->gamma ) ) this->positiveMainMomentumFlux( u[ center ], velocity_x_center, pressure_center )
* ( - this->positiveMainMomentumFlux( u[ west ], velocity_x_west , pressure_west )
2 * ( this->gamma - 1 ) * velocity_x_center * velocity_x_center - this->positiveMainMomentumFlux( u[ center ], velocity_x_center, pressure_center )
+ ( velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) ) + this->positiveMainMomentumFlux( u[ east ], velocity_x_east , pressure_east )
* ( 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 ] ) )
); );
} }
......
examples/inviscid-flow/RiemannProblemInitialCondition.h
+ 1
− 1
View file @ 8286beec
...@@ -949,7 +949,7 @@ class RiemannProblemInitialCondition ...@@ -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 preNEUVelocityX, double preNEUVelocityY,double preNEUVelocityZ,
0.0, 0.0, 0.0, //double preSEUVelocityX, double preSEUVelocityY,double preSEUVelocityZ, 0.0, 0.0, 0.0, //double preSEUVelocityX, double preSEUVelocityY,double preSEUVelocityZ,
0.0, 0.0, 0.0, //double preNEDVelocityX, double preNEDVelocityY,double preNEDVelocityZ, 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") if(initial == prefix + "1D_4")
predefinedInitialCondition( 1.666, 0.4, 0.0, 0.0, // double preGamma, double preDiscX, double preDiscY, double preDiscZ, predefinedInitialCondition( 1.666, 0.4, 0.0, 0.0, // double preGamma, double preDiscX, double preDiscY, double preDiscZ,
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment