From 8286beeca8d71ba112e24ffaed2c4232a9dc94ed Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jan=20Sch=C3=A4fer?= <schafjan@fjfi.cvut.cz>
Date: Sat, 12 May 2018 17:37:47 +0200
Subject: [PATCH] plne funkcni 1D, odzkouseno
---
.../RiemannProblemInitialCondition.h | 2 +-
examples/inviscid-flow-sw/UpwindEnergy.h | 92 +++++--------------
.../inviscid-flow-sw/UpwindMomentumBase.h | 16 ++--
examples/inviscid-flow-sw/UpwindMomentumX.h | 22 +----
examples/inviscid-flow-sw/UpwindMomentumX.h~ | 22 +----
.../RiemannProblemInitialCondition.h | 2 +-
6 files changed, 43 insertions(+), 113 deletions(-)
diff --git a/examples/inviscid-flow-sw/RiemannProblemInitialCondition.h b/examples/inviscid-flow-sw/RiemannProblemInitialCondition.h
index bfa458a8bb..85fbc8af01 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 7ecc6449c2..29d74e32e9 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 c9bb5745fc..29d4612c1d 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 662e0707f4..c9cafec441 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 4b941f23ae..7f72b31dae 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 a036747f25..b1819b1eb8 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,
--
GitLab