nove resen upvind pomoci stegger warming zatim nekompletni. (3008a6e7) · Commits · TNL / tnl-dev

examples/inviscid-flow-sw/RiemannProblemInitialCondition.h

+2 −2

Original line number	Diff line number	Diff line
		@@ -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,

examples/inviscid-flow-sw/Upwind.h

+9 −2

Original line number	Diff line number	Diff line
		@@ -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;
		}

		@@ -79,6 +77,15 @@ class Upwind
		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 )
		{
		this->continuityOperatorPointer->setPressure( pressure );

examples/inviscid-flow-sw/UpwindContinuity.h

+45 −13

Original line number	Diff line number	Diff line
		@@ -65,6 +65,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:

		@@ -134,10 +162,14 @@ 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 ] ) )
		// - 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 ] ) )
		);

		}

examples/inviscid-flow-sw/UpwindEnergy.h

+60 −10

Original line number	Diff line number	Diff line
		@@ -69,6 +69,54 @@ 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 )
		{
		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:

		RealType tau;
		@@ -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 * 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 * pressure_east / u[ east ] )
		/ 2
		)
		)*/
		);

		}

examples/inviscid-flow-sw/UpwindMomentumBase.h

+56 −0

Original line number	Diff line number	Diff line
		@@ -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;