From 3008a6e72ef2e3881b1f09b6358f00c156414222 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jan=20Sch=C3=A4fer?= <schafjan@fjfi.cvut.cz>
Date: Wed, 9 May 2018 20:05:57 +0200
Subject: [PATCH] nove resen upvind pomoci stegger warming zatim nekompletni.
---
.../RiemannProblemInitialCondition.h | 4 +-
examples/inviscid-flow-sw/Upwind.h | 11 ++-
examples/inviscid-flow-sw/UpwindContinuity.h | 58 +++++++++++----
examples/inviscid-flow-sw/UpwindEnergy.h | 70 ++++++++++++++++---
.../inviscid-flow-sw/UpwindMomentumBase.h | 56 +++++++++++++++
examples/inviscid-flow-sw/UpwindMomentumX.h | 12 ++--
examples/inviscid-flow-sw/UpwindMomentumX.h~ | 60 ++++++++--------
examples/inviscid-flow-sw/eulerProblem_impl.h | 23 ++++--
8 files changed, 224 insertions(+), 70 deletions(-)
diff --git a/examples/inviscid-flow-sw/RiemannProblemInitialCondition.h b/examples/inviscid-flow-sw/RiemannProblemInitialCondition.h
index a036747f25..bfa458a8bb 100644
--- a/examples/inviscid-flow-sw/RiemannProblemInitialCondition.h
+++ b/examples/inviscid-flow-sw/RiemannProblemInitialCondition.h
@@ -990,11 +990,11 @@ class RiemannProblemInitialCondition
0.0, 0.0, 0.0, //double preNWUVelocityX, double preNWUVelocityY,double preNWUVelocityZ,
0.0, 0.0, 0.0, //double preSWUVelocityX, double preSWUVelocityY,double preSWUVelocityZ,
0.0, 0.0, 0.0, //double preNWDVelocityX, double preNWDVelocityY,double preNWDVelocityZ,
- 1.0, 0.0, 0.0, //double preSWDVelocityX, double preSWDVelocityY,double preSWDVelocityZ,
+ 0.1, 0.0, 0.0, //double preSWDVelocityX, double preSWDVelocityY,double preSWDVelocityZ,
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,
- 1.0, 0.0, 0.0 //double preSEDVelocityX, double preSEDVelocityY,double preSEDVelocityZ
+ 0.1, 0.0, 0.0 //double preSEDVelocityX, double preSEDVelocityY,double preSEDVelocityZ
);
if(initial == prefix + "1D_Noh")
predefinedInitialCondition( 1.4, 0.5, 0.0, 0.0, // double preGamma, double preDiscX, double preDiscY, double preDiscZ,
diff --git a/examples/inviscid-flow-sw/Upwind.h b/examples/inviscid-flow-sw/Upwind.h
index 3e2140da8f..38776d80a0 100644
--- a/examples/inviscid-flow-sw/Upwind.h
+++ b/examples/inviscid-flow-sw/Upwind.h
@@ -56,7 +56,6 @@ class Upwind
static void configSetup( Config::ConfigDescription& config,
const String& prefix = "" )
{
- config.addEntry< double >( prefix + "numerical-viscosity", "Value of artificial (numerical) viscosity in the Lax-Fridrichs scheme", 1.0 );
}
Upwind()
@@ -66,7 +65,6 @@ class Upwind
const Config::ParameterContainer& parameters,
const String& prefix = "" )
{
- this->artificialViscosity = parameters.getParameter< double >( prefix + "numerical-viscosity" );
return true;
}
@@ -78,6 +76,15 @@ class Upwind
this->momentumZOperatorPointer->setTau( tau );
this->energyOperatorPointer->setTau( tau );
}
+
+ void setGamma( const RealType& gamma )
+ {
+ this->continuityOperatorPointer->setGamma( gamma );
+ this->momentumXOperatorPointer->setGamma( gamma );
+ this->momentumYOperatorPointer->setGamma( gamma );
+ this->momentumZOperatorPointer->setGamma( gamma );
+ this->energyOperatorPointer->setGamma( gamma );
+ }
void setPressure( const MeshFunctionPointer& pressure )
{
diff --git a/examples/inviscid-flow-sw/UpwindContinuity.h b/examples/inviscid-flow-sw/UpwindContinuity.h
index 5184163ae8..4439d83138 100644
--- a/examples/inviscid-flow-sw/UpwindContinuity.h
+++ b/examples/inviscid-flow-sw/UpwindContinuity.h
@@ -64,6 +64,34 @@ class UpwindContinuityBase
{
this->velocity = velocity;
};
+
+ const RealType& positiveDensityFlux( 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 0.0;
+ else if ( machNumber <= 0.0 )
+ return density * speedOfSound / ( 2 * this->gamma ) * ( machNumber + 1.0 );
+ else if ( machNumber <= 1.0 )
+ 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 )
+ {
+ const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
+ const RealType& machNumber = velocity / speedOfSound;
+ if ( machNumber <= -1.0 )
+ return density * velocity;
+ else if ( machNumber <= 0.0 )
+ return density * speedOfSound / ( 2 * this->gamma ) * ( ( 2.0 * this->gamma - 1.0 ) * machNumber - 1.0 );
+ else if ( machNumber <= 1.0 )
+ return density * speedOfSound / ( 2 * this->gamma ) * ( machNumber - 1.0 );
+ else
+ return 0.0;
+ }
protected:
@@ -133,11 +161,15 @@ class UpwindContinuity< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, In
const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
const RealType& velocity_x_east = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ];
- return -hxInverse * (
+ 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 ] ) )
// - 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 ] ) )
+// - 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 ] ) )
);
}
@@ -217,14 +249,14 @@ class UpwindContinuity< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, In
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 ] * ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) ) / ( 2 * this->gamma )
- + u[ east ] * ( velocity_x_east - std::sqrt( this->gamma * pressure_east / u[ east ] ) ) / ( 2 * this->gamma )
+ - 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 ] ) )
)
-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 ] * ( velocity_y_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) ) / ( 2 * this->gamma )
- + u[ north ] * ( velocity_y_north - std::sqrt( this->gamma * pressure_north / u[ north ] ) ) / ( 2 * this->gamma )
+ - 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 ] ) )
);
}
@@ -313,20 +345,20 @@ class UpwindContinuity< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, In
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 ] * ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) ) / ( 2 * this->gamma )
- + u[ east ] * ( velocity_x_east - std::sqrt( this->gamma * pressure_east / u[ east ] ) ) / ( 2 * this->gamma )
+ - 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 ] ) )
)
-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 ] * ( velocity_y_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) ) / ( 2 * this->gamma )
- + u[ north ] * ( velocity_y_north - std::sqrt( this->gamma * pressure_north / u[ north ] ) ) / ( 2 * this->gamma )
+ - 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 ] ) )
)
-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 ] * ( velocity_z_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) ) / ( 2 * this->gamma )
- + u[ up ] * ( velocity_z_up - std::sqrt( this->gamma * pressure_up / u[ up ] ) ) / ( 2 * this->gamma )
+ - 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 ] ) )
);
}
diff --git a/examples/inviscid-flow-sw/UpwindEnergy.h b/examples/inviscid-flow-sw/UpwindEnergy.h
index 494fdbbe56..8029f3dc69 100644
--- a/examples/inviscid-flow-sw/UpwindEnergy.h
+++ b/examples/inviscid-flow-sw/UpwindEnergy.h
@@ -67,6 +67,54 @@ 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 )
+ {
+ 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 )
+ * (
+ ( 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 * 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 );
+ }
+
+ const RealType& negativeEnergyFlux( 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 ( 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 * 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 * speedOfSound * ( machNumber - 1.0 ) / ( this->gamma - 1.0 ) )
+ );
+ else
+ return 0.0;
}
protected:
@@ -151,36 +199,38 @@ class UpwindEnergy< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index
)
- ( u[ west ] / ( 2 * this->gamma ) )
* (
- ( 2 * this->gamma - 1 ) * velocity_x_west * velocity_x_west
+ ( 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_center / u[ west ] ) )
+ * ( velocity_x_west + std::sqrt( this->gamma * pressure_west / u[ west ] ) )
* ( this->gamma * pressure_west / u[ west ] )
)
- - u[ center ] / ( 2 * this->gamma )
+/* - 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 )
- * ( velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
+ / ( 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_center / u[ east ] ) )
- * ( velocity_x_east - std::sqrt( this->gamma * pressure_center / 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 )
- * ( velocity_x_east + std::sqrt( this->gamma * pressure_east / u[ east ] ) )
+ / ( this->gamma - 1 )
+ * ( this->gamma * pressure_east / u[ east ] )
/ 2
- )
+ )*/
);
}
diff --git a/examples/inviscid-flow-sw/UpwindMomentumBase.h b/examples/inviscid-flow-sw/UpwindMomentumBase.h
index 818d6cc56d..8bbc9015a9 100644
--- a/examples/inviscid-flow-sw/UpwindMomentumBase.h
+++ b/examples/inviscid-flow-sw/UpwindMomentumBase.h
@@ -52,6 +52,62 @@ class UpwindMomentumBase
this->pressure = pressure;
};
+ const RealType& positiveMainMomentumFlux( 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 0;
+ else if ( machNumber <= 0.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 );
+ 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;
+ 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 );
+ else if ( machNumber <= 1.0 )
+ return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( ( machNumber - 1.0 ) * machNumber + ( - machNumber + 1.0 );
+ else
+ return density * velocity * velocity * pressure;
+ }
+
+ const RealType& positiveOtherMomentumFlux( 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 0;
+ else if ( machNumber <= 0.0 )
+ return density * speedOfSound / ( 2 * this->gamma ) * ( machNumber + 1.0 );
+ else if ( machNumber <= 1.0 )
+ return density * speedOfSound / ( 2 * this->gamma ) * ( ( 2.0 * this->gamma - 1.0 ) * machNumber + 1.0 );
+ else
+ return density * velocity;
+ }
+
+ const RealType& negativeOtherMomentumFlux( 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;
+ else if ( machNumber <= 0.0 )
+ return density * speedOfSound / ( 2 * this->gamma ) * ( ( 2.0 * this->gamma - 1.0 ) * machNumber - 1.0 );
+ else if ( machNumber <= 1.0 )
+ return density * speedOfSound / ( 2 * this->gamma ) * ( machNumber - 1.0 );
+ else
+ return density * velocity;
+ }
+
protected:
RealType tau;
diff --git a/examples/inviscid-flow-sw/UpwindMomentumX.h b/examples/inviscid-flow-sw/UpwindMomentumX.h
index 347306fa8a..662e0707f4 100644
--- a/examples/inviscid-flow-sw/UpwindMomentumX.h
+++ b/examples/inviscid-flow-sw/UpwindMomentumX.h
@@ -90,12 +90,12 @@ class UpwindMomentumX< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Ind
+ ( velocity_x_west + std::sqrt( this->gamma * pressure_west / u[ west ] ) )
* ( velocity_x_west + std::sqrt( this->gamma * pressure_west / u[ west ] ) )
)
- - u[ center ]
- * ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) ) / ( 2 * this->gamma )
- * ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) ) / ( 2 * this->gamma )
- + u[ east ]
- * ( velocity_x_east - std::sqrt( this->gamma * pressure_east / u[ east ] ) ) / ( 2 * this->gamma )
- * ( velocity_x_east - std::sqrt( this->gamma * pressure_east / u[ east ] ) ) / ( 2 * this->gamma )
+/* - ( 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 ] ) )*/
);
}
diff --git a/examples/inviscid-flow-sw/UpwindMomentumX.h~ b/examples/inviscid-flow-sw/UpwindMomentumX.h~
index 9b20007e1d..4b941f23ae 100644
--- a/examples/inviscid-flow-sw/UpwindMomentumX.h~
+++ b/examples/inviscid-flow-sw/UpwindMomentumX.h~
@@ -58,7 +58,7 @@ class UpwindMomentumX< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Ind
template< typename MeshFunction, typename MeshEntity >
__cuda_callable__
- Real operator()( const MeshFunction& rho_u,
+ Real operator()( const MeshFunction& u,
const MeshEntity& entity,
const RealType& time = 0.0 ) const
{
@@ -70,34 +70,32 @@ 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& rho_center = this->rho.template getData< DeviceType >()[ center ];
- const RealType& rho_west = this->rho.template getData< DeviceType >()[ west ];
- const RealType& rho_east = this->rho.template getData< DeviceType >()[ east ];
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 * (
( u[ center ] / ( 2 * this->gamma ) )
* (
- ( 2 * this->gamma - 1 ) * velocity_x_center * velocity_x_center
+ 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
+ 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 ]
- * ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) ) / ( 2 * this->gamma )
- * ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) ) / ( 2 * this->gamma )
- + u[ east ]
- * ( velocity_x_east - std::sqrt( this->gamma * pressure_east / u[ east ] ) ) / ( 2 * this->gamma )
- * ( velocity_x_east - std::sqrt( this->gamma * pressure_east / u[ east ] ) ) / ( 2 * this->gamma )
+ - ( 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 ] ) )
);
}
@@ -151,7 +149,7 @@ class UpwindMomentumX< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Ind
template< typename MeshFunction, typename MeshEntity >
__cuda_callable__
- Real operator()( const MeshFunction& rho_u,
+ Real operator()( const MeshFunction& u,
const MeshEntity& entity,
const RealType& time = 0.0 ) const
{
@@ -189,12 +187,12 @@ class UpwindMomentumX< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Ind
return -hxInverse * (
u[ center ] / ( 2 * this->gamma )
- * ( ( 2 * this->gamma - 1 ) * velocity_x_center * velocity_x_center
+ * ( 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
+ * ( 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 ] ) )
)
@@ -207,16 +205,16 @@ class UpwindMomentumX< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Ind
)
-hyInverse * (
u[ center ] / ( 2 * this->gamma )
- * ( ( 2 * this->gamma - 1 ) * this->velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
+ * ( ( 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 ) * this->velocity_x_south + std::sqrt( this->gamma * pressure_south / u[ south ] ) )
+ * ( ( 2 * this->gamma - 1 ) * velocity_x_south + std::sqrt( this->gamma * pressure_south / u[ south ] ) )
* velocity_y_south
- u[ center ] / ( 2 * this->gamma )
- * ( this->velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
+ * ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
* velocity_y_center
+ u[ north ] / ( 2 * this->gamma )
- * ( this->velocity_x_north - std::sqrt( this->gamma * pressure_north / u[ north ] ) )
+ * ( velocity_x_north - std::sqrt( this->gamma * pressure_north / u[ north ] ) )
* velocity_y_north
);
@@ -273,7 +271,7 @@ class UpwindMomentumX< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real, Inde
template< typename MeshFunction, typename MeshEntity >
__cuda_callable__
- Real operator()( const MeshFunction& rho_u,
+ Real operator()( const MeshFunction& u,
const MeshEntity& entity,
const RealType& time = 0.0 ) const
{
@@ -313,7 +311,7 @@ class UpwindMomentumX< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real, Inde
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 >()[ souh ];
+ 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 ];
@@ -327,12 +325,12 @@ class UpwindMomentumX< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real, Inde
return -hxInverse * (
u[ center ] / ( 2 * this->gamma )
- * ( ( 2 * this->gamma - 1 ) * velocity_x_center * velocity_x_center
+ * ( 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
+ * ( 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 ] ) )
)
@@ -345,30 +343,30 @@ class UpwindMomentumX< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real, Inde
)
-hyInverse * (
u[ center ] / ( 2 * this->gamma )
- * ( ( 2 * this->gamma - 1 ) * this->velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
+ * ( 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 ) * this->velocity_x_south + std::sqrt( this->gamma * pressure_south / u[ south ] ) )
+ * ( 2 * ( this->gamma - 1 ) * velocity_x_south + std::sqrt( this->gamma * pressure_south / u[ south ] ) )
* velocity_y_south
- u[ center ] / ( 2 * this->gamma )
- * ( this->velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
+ * ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
* velocity_y_center
+ u[ north ] / ( 2 * this->gamma )
- * ( this->velocity_x_north - std::sqrt( this->gamma * pressure_north / u[ north ] ) )
+ * ( velocity_x_north - std::sqrt( this->gamma * pressure_north / u[ north ] ) )
* velocity_y_north
)
-hzInverse * (
u[ center ] / ( 2 * this->gamma )
- * ( ( 2 * this->gamma - 1 ) * this->velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
+ * ( 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 ) * this->velocity_x_down + std::sqrt( this->gamma * pressure_down / u[ down ] ) )
+ * ( 2 * ( this->gamma - 1 ) * velocity_x_down + std::sqrt( this->gamma * pressure_down / u[ down ] ) )
* velocity_z_south
- u[ center ] / ( 2 * this->gamma )
- * ( this->velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
+ * ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )
* velocity_z_center
+ u[ up ] / ( 2 * this->gamma )
- * ( this->velocity_x_up - std::sqrt( this->gamma * pressure_up / u[ up ] ) )
+ * ( velocity_x_up - std::sqrt( this->gamma * pressure_up / u[ up ] ) )
* velocity_z_up
);
}
diff --git a/examples/inviscid-flow-sw/eulerProblem_impl.h b/examples/inviscid-flow-sw/eulerProblem_impl.h
index 2e2cda9a25..e1d5ea9fb1 100644
--- a/examples/inviscid-flow-sw/eulerProblem_impl.h
+++ b/examples/inviscid-flow-sw/eulerProblem_impl.h
@@ -199,9 +199,13 @@ makeSnapshot( const RealType& time,
if( ! this->pressure->save( fileName.getFileName() ) )
return false;
-// fileName.setFileNameBase( "energy-" );
-// if( ! this->conservativeVariables->getEnergy()->save( fileName.getFileName() ) )
-// return false;
+ fileName.setFileNameBase( "energy-" );
+ if( ! this->conservativeVariables->getEnergy()->save( fileName.getFileName() ) )
+ return false;
+
+ fileName.setFileNameBase( "momentum-" );
+ if( ! this->conservativeVariables->getMomentum()->save( fileName.getFileName() ) )
+ return false;
return true;
}
@@ -228,7 +232,9 @@ getExplicitUpdate( const RealType& time,
this->conservativeVariablesRHS->bind( mesh, _fu );
this->velocity->setMesh( mesh );
this->pressure->setMesh( mesh );
-
+
+// this->pressure->write( "pressure1", "gnuplot" );
+// getchar();
/****
* Resolve the physical variables
*/
@@ -248,7 +254,10 @@ getExplicitUpdate( const RealType& time,
this->inviscidOperatorsPointer->setTau( tau );
this->inviscidOperatorsPointer->setVelocity( this->velocity );
this->inviscidOperatorsPointer->setPressure( this->pressure );
+ this->inviscidOperatorsPointer->setGamma( this->gamma );
+// this->pressure->write( "pressure2", "gnuplot" );
+// getchar();
/****
* Continuity equation
*/
@@ -304,8 +313,10 @@ getExplicitUpdate( const RealType& time,
explicitUpdaterEnergy.template update< typename Mesh::Cell >( time, tau, mesh,
this->conservativeVariables->getDensity(), // uRhoVelocityX,
this->conservativeVariablesRHS->getEnergy() ); //, fuRhoVelocityX );
-
- /*this->conservativeVariablesRHS->getDensity()->write( "density", "gnuplot" );
+
+/* this->pressure->write( "pressure3", "gnuplot" );
+ getchar();
+ this->conservativeVariablesRHS->getDensity()->write( "density", "gnuplot" );
this->conservativeVariablesRHS->getEnergy()->write( "energy", "gnuplot" );
this->conservativeVariablesRHS->getMomentum()->write( "momentum", "gnuplot", 0.05 );
getchar();*/
--
GitLab