Adding artificial viscosity to the Lax-Fridrichs scheme for inviscid flow. (2e18d455) · Commits · TNL / tnl-dev

examples/inviscid-flow/LaxFridrichs.h

+103 −14

Original line number	Diff line number	Diff line
		@@ -13,16 +13,13 @@

		#include <TNL/Containers/Vector.h>
		#include <TNL/Meshes/Grid.h>
		#include <TNL/Functions/VectorField.h>

		#include "LaxFridrichsContinuity.h"
		#include "LaxFridrichsEnergy.h"
		#include "LaxFridrichsMomentumX.h"
		#include "LaxFridrichsMomentumY.h"
		#include "LaxFridrichsMomentumZ.h"
		/*#include "EulerPressureGetter.h"
		#include "EulerVelXGetter.h"
		#include "EulerVelYGetter.h"
		#include "EulerVelGetter.h"*/

		namespace TNL {

		@@ -32,21 +29,113 @@ template< typename Mesh,
		class LaxFridrichs
		{
		public:
		typedef Mesh MeshType;
		typedef Real RealType;
		typedef typename Mesh::DeviceType DeviceType;
		typedef Index IndexType;
		typedef Functions::MeshFunction< Mesh > MeshFunctionType;
		static const int Dimensions = Mesh::getMeshDimensions();
		typedef Functions::VectorField< Dimensions, MeshFunctionType > VectorFieldType;

		typedef LaxFridrichsContinuity< Mesh, Real, Index > Continuity;
		typedef LaxFridrichsMomentumX< Mesh, Real, Index > MomentumX;
		typedef LaxFridrichsMomentumY< Mesh, Real, Index > MomentumY;
		typedef LaxFridrichsMomentumZ< Mesh, Real, Index > MomentumZ;
		typedef LaxFridrichsEnergy< Mesh, Real, Index > Energy;
		/*typedef EulerVelXGetter< Mesh, Real, Index > VelocityX;
		typedef EulerVelYGetter< Mesh, Real, Index > VelocityY;
		typedef EulerVelGetter< Mesh, Real, Index > Velocity;
		typedef EulerPressureGetter< Mesh, Real, Index > Pressure;*/
		typedef LaxFridrichsContinuity< Mesh, Real, Index > ContinuityOperatorType;
		typedef LaxFridrichsMomentumX< Mesh, Real, Index > MomentumXOperatorType;
		typedef LaxFridrichsMomentumY< Mesh, Real, Index > MomentumYOperatorType;
		typedef LaxFridrichsMomentumZ< Mesh, Real, Index > MomentumZOperatorType;
		typedef LaxFridrichsEnergy< Mesh, Real, Index > EnergyOperatorType;

		typedef SharedPointer< MeshFunctionType > MeshFunctionPointer;
		typedef SharedPointer< VectorFieldType > VectorFieldPointer;
		typedef SharedPointer< MeshType > MeshPointer;

		typedef SharedPointer< ContinuityOperatorType > ContinuityOperatorPointer;
		typedef SharedPointer< MomentumXOperatorType > MomentumXOperatorPointer;
		typedef SharedPointer< MomentumYOperatorType > MomentumYOperatorPointer;
		typedef SharedPointer< MomentumZOperatorType > MomentumZOperatorPointer;
		typedef SharedPointer< EnergyOperatorType > EnergyOperatorPointer;

		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 );
		}

		LaxFridrichs()
		: artificialViscosity( 1.0 ) {}

		bool setup( const MeshPointer& meshPointer,
		const Config::ParameterContainer& parameters,
		const String& prefix = "" )
		{
		this->artificialViscosity = parameters.getParameter< double >( prefix + "numerical-viscosity" );
		this->continuityOperatorPointer->setArtificialViscosity( artificialViscosity );
		this->momentumXOperatorPointer->setArtificialViscosity( artificialViscosity );
		this->momentumYOperatorPointer->setArtificialViscosity( artificialViscosity );
		this->momentumZOperatorPointer->setArtificialViscosity( artificialViscosity );
		this->energyOperatorPointer->setArtificialViscosity( artificialViscosity );

		return true;
		}

		void setTau( const RealType& tau )
		{
		this->continuityOperatorPointer->setTau( tau );
		this->momentumXOperatorPointer->setTau( tau );
		this->momentumYOperatorPointer->setTau( tau );
		this->momentumZOperatorPointer->setTau( tau );
		this->energyOperatorPointer->setTau( tau );
		}

		void setPressure( const MeshFunctionPointer& pressure )
		{
		this->momentumXOperatorPointer->setPressure( pressure );
		this->momentumYOperatorPointer->setPressure( pressure );
		this->momentumZOperatorPointer->setPressure( pressure );
		this->energyOperatorPointer->setPressure( pressure );
		}

		void setVelocity( const VectorFieldPointer& velocity )
		{
		this->continuityOperatorPointer->setVelocity( velocity );
		this->momentumXOperatorPointer->setVelocity( velocity );
		this->momentumYOperatorPointer->setVelocity( velocity );
		this->momentumZOperatorPointer->setVelocity( velocity );
		this->energyOperatorPointer->setVelocity( velocity );
		}

		const ContinuityOperatorPointer& getContinuityOperator() const
		{
		return this->continuityOperatorPointer;
		}

		const MomentumXOperatorPointer& getMomentumXOperator() const
		{
		return this->momentumXOperatorPointer;
		}

		const MomentumYOperatorPointer& getMomentumYOperator() const
		{
		return this->momentumYOperatorPointer;
		}

		const MomentumZOperatorPointer& getMomentumZOperator() const
		{
		return this->momentumZOperatorPointer;
		}

		const EnergyOperatorPointer& getEnergyOperator() const
		{
		return this->energyOperatorPointer;
		}

		protected:

		ContinuityOperatorPointer continuityOperatorPointer;
		MomentumXOperatorPointer momentumXOperatorPointer;
		MomentumYOperatorPointer momentumYOperatorPointer;
		MomentumZOperatorPointer momentumZOperatorPointer;
		EnergyOperatorPointer energyOperatorPointer;

		RealType artificialViscosity;
		};

		} //namespace TNL

examples/inviscid-flow/LaxFridrichsContinuity.h

+16 −5

Original line number	Diff line number	Diff line
		@@ -36,6 +36,9 @@ class LaxFridrichsContinuityBase
		typedef Functions::VectorField< Dimensions, MeshFunctionType > VelocityFieldType;
		typedef SharedPointer< VelocityFieldType > VelocityFieldPointer;

		LaxFridrichsContinuityBase()
		: artificialViscosity( 1.0 ){};

		static String getType()
		{
		return String( "LaxFridrichsContinuity< " ) +
		@@ -54,12 +57,19 @@ class LaxFridrichsContinuityBase
		this->velocity = velocity;
		};

		void setArtificialViscosity( const RealType& artificialViscosity )
		{
		this->artificialViscosity = artificialViscosity;
		}


		protected:

		RealType tau;

		VelocityFieldPointer velocity;

		RealType artificialViscosity;
		};


		@@ -91,7 +101,7 @@ class LaxFridrichsContinuity< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Re
		using typename BaseType::MeshFunctionType;
		using typename BaseType::VelocityFieldType;
		using typename BaseType::VelocityFieldPointer;
		using typename BaseType::Dimensions;
		using BaseType::Dimensions;

		template< typename MeshFunction, typename MeshEntity >
		__cuda_callable__
		@@ -109,7 +119,7 @@ class LaxFridrichsContinuity< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Re
		const IndexType& west = neighbourEntities.template getEntityIndex< -1 >();
		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 1.0 / 2.0 * this->tau * ( u[ west ] - 2.0 * u[ center ] + u[ east ] )
		return 1.0 / 2.0 * this->tau * this->artificialViscosity * ( u[ west ] - 2.0 * u[ center ] + u[ east ] )
		- 0.5 * ( u[ west ] * velocity_x_west - u[ east ] * velocity_x_east ) * hxInverse;
		}

		@@ -175,7 +185,7 @@ class LaxFridrichsContinuity< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Re
		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 >()[ south ];

		return 1.0 / 4.0 * this->tau * ( u[ west ] + u[ east ] + u[ south ] + u[ north ] - 4.0 * u[ center ] )
		return 1.0 / 4.0 * this->tau * this->artificialViscosity * ( u[ west ] + u[ east ] + u[ south ] + u[ north ] - 4.0 * u[ center ] )
		- 0.5 * ( ( u[ west ] * velocity_x_west - u[ east ] * velocity_x_east ) * hxInverse
		+ ( u[ north ] * velocity_y_north - u[ south ] * velocity_y_south ) * hyInverse );
		}
		@@ -213,11 +223,11 @@ class LaxFridrichsContinuity< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Re
		using typename BaseType::RealType;
		using typename BaseType::IndexType;
		using typename BaseType::DeviceType;
		using typename BaseType::Dimensions;
		using typename BaseType::CoordinatesType;
		using typename BaseType::MeshFunctionType;
		using typename BaseType::VelocityFieldType;
		using typename BaseType::VelocityFieldPointer;
		using BaseType::Dimensions;

		template< typename MeshFunction, typename MeshEntity >
		__cuda_callable__
		@@ -248,7 +258,8 @@ class LaxFridrichsContinuity< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Re
		const RealType& velocity_z_up = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ up ];
		const RealType& velocity_z_down = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ down ];

		return 1.0 / 6.0 * this->tau * ( u[ west ] + u[ east ] + u[ south ] + u[ north ] + u[ up ] + u[ down ]- 6.0 * u[ center ] )
		return 1.0 / 6.0 * this->tau * this->artificialViscosity *
		( u[ west ] + u[ east ] + u[ south ] + u[ north ] + u[ up ] + u[ down ]- 6.0 * u[ center ] )
		- 0.5 * ( ( u[ west ] * velocity_x_west - u[ east ] * velocity_x_east ) * hxInverse
		+ ( u[ north ] * velocity_y_north - u[ south ] * velocity_y_south ) * hyInverse
		+ ( u[ up ] * velocity_z_up - u[ down ] * velocity_z_down ) * hzInverse );

examples/inviscid-flow/LaxFridrichsEnergy.h

+14 −3

Original line number	Diff line number	Diff line
		@@ -33,6 +33,9 @@ class LaxFridrichsEnergyBase
		typedef SharedPointer< MeshFunctionType > MeshFunctionPointer;
		typedef SharedPointer< VelocityFieldType > VelocityFieldPointer;

		LaxFridrichsEnergyBase()
		: artificialViscosity( 1.0 ){};

		static String getType()
		{
		return String( "LaxFridrichsEnergy< " ) +
		@@ -56,6 +59,11 @@ class LaxFridrichsEnergyBase
		this->pressure = pressure;
		};

		void setArtificialViscosity( const RealType& artificialViscosity )
		{
		this->artificialViscosity = artificialViscosity;
		}

		protected:

		RealType tau;
		@@ -63,6 +71,8 @@ class LaxFridrichsEnergyBase
		VelocityFieldPointer velocity;

		MeshFunctionPointer pressure;

		RealType artificialViscosity;
		};

		template< typename Mesh,
		@@ -113,7 +123,7 @@ class LaxFridrichsEnergy< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real,
		const RealType& pressure_east = this->pressure.template getData< DeviceType >()[ east ];
		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 1.0 / 2.0 * this->tau * ( e[ west ] - 2.0 * e[ center ] + e[ east ] )
		return 1.0 / 2.0 * this->tau * this->artificialViscosity * ( e[ west ] - 2.0 * e[ center ] + e[ east ] )
		- 0.5 * ( ( e[ west ] + pressure_west ) * velocity_x_west
		- ( e[ east ] + pressure_east ) * velocity_x_east ) * hxInverse;

		@@ -186,7 +196,7 @@ class LaxFridrichsEnergy< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real,
		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 >()[ south ];

		return 1.0 / 4.0 * this->tau * ( e[ west ] + e[ east ] + e[ south ] + e[ north ] - 4.0 * e[ center ] )
		return 1.0 / 4.0 * this->tau * this->artificialViscosity * ( e[ west ] + e[ east ] + e[ south ] + e[ north ] - 4.0 * e[ center ] )
		- 0.5 * ( ( ( ( e[ west ] + pressure_west ) * velocity_x_west )
		-( ( e[ east ] + pressure_east ) * velocity_x_east ) ) * hxInverse
		+ ( ( ( e[ north ] + pressure_north ) * velocity_y_north )
		@@ -268,7 +278,8 @@ class LaxFridrichsEnergy< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real,
		const RealType& velocity_z_up = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ up ];
		const RealType& velocity_z_down = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ down ];

		return 1.0 / 6.0 * this->tau * ( e[ west ] + e[ east ] + e[ south ] + e[ north ] + e[ up ] + e[ down ] - 6.0 * e[ center ] )
		return 1.0 / 6.0 * this->tau * this->artificialViscosity *
		( e[ west ] + e[ east ] + e[ south ] + e[ north ] + e[ up ] + e[ down ] - 6.0 * e[ center ] )
		- 0.5 * ( ( ( ( e[ west ] + pressure_west ) * velocity_x_west )
		-( ( e[ east ] + pressure_east ) * velocity_x_east ) ) * hxInverse
		+ ( ( ( e[ north ] + pressure_north ) * velocity_y_north )

examples/inviscid-flow/LaxFridrichsMomentumBase.h

+9 −0

Original line number	Diff line number	Diff line
		@@ -31,6 +31,8 @@ class LaxFridrichsMomentumBase
		typedef SharedPointer< MeshFunctionType > MeshFunctionPointer;
		typedef SharedPointer< VelocityFieldType > VelocityFieldPointer;

		LaxFridrichsMomentumBase()
		: artificialViscosity( 1.0 ){};

		void setTau(const Real& tau)
		{
		@@ -47,6 +49,11 @@ class LaxFridrichsMomentumBase
		this->pressure = pressure;
		};

		void setArtificialViscosity( const RealType& artificialViscosity )
		{
		this->artificialViscosity = artificialViscosity;
		}

		protected:

		RealType tau;
		@@ -54,6 +61,8 @@ class LaxFridrichsMomentumBase
		VelocityFieldPointer velocity;

		MeshFunctionPointer pressure;

		RealType artificialViscosity;
		};

		} //namespace TNL

examples/inviscid-flow/LaxFridrichsMomentumX.h

+4 −3

Original line number	Diff line number	Diff line
		@@ -75,7 +75,7 @@ class LaxFridrichsMomentumX< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Rea
		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 1.0 / 2.0 * this->tau * ( rho_u[ west ] + rho_u[ east ] - 2.0 * rho_u[ center ] )
		return 1.0 / 2.0 * this->tau * this->artificialViscosity * ( rho_u[ west ] + rho_u[ east ] - 2.0 * rho_u[ center ] )
		- 0.5 * ( ( rho_u[ west ] * velocity_x_west + pressure_west )
		-( rho_u[ east ] * velocity_x_east + pressure_east ) ) * hxInverse;
		}
		@@ -154,7 +154,7 @@ class LaxFridrichsMomentumX< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Rea
		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 >()[ south ];

		return 1.0 / 4.0 * this->tau * ( rho_u[ west ] + rho_u[ east ] + rho_u[ south ] + rho_u[ north ] - 4.0 * rho_u[ center ] )
		return 1.0 / 4.0 * this->tau * this->artificialViscosity * ( rho_u[ west ] + rho_u[ east ] + rho_u[ south ] + rho_u[ north ] - 4.0 * rho_u[ center ] )
		- 0.5 * ( ( ( rho_u[ west ] * velocity_x_west + pressure_west )
		- ( rho_u[ east ] * velocity_x_east + pressure_east ) ) * hxInverse
		+ ( ( rho_u[ north ] * velocity_y_north )
		@@ -243,7 +243,8 @@ class LaxFridrichsMomentumX< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real
		const RealType& velocity_y_south = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ];
		const RealType& velocity_z_up = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ up ];
		const RealType& velocity_z_down = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ down ];
		return 1.0 / 6.0 * this->tau * ( rho_u[ west ] + rho_u[ east ] + rho_u[ south ] + rho_u[ north ] + rho_u[ up ] + rho_u[ down ] - 6.0 * rho_u[ center ] )
		return 1.0 / 6.0 * this->tau * this->artificialViscosity *
		( rho_u[ west ] + rho_u[ east ] + rho_u[ south ] + rho_u[ north ] + rho_u[ up ] + rho_u[ down ] - 6.0 * rho_u[ center ] )
		- 0.5 * ( ( ( rho_u[ west ] * velocity_x_west + pressure_west )
		- ( rho_u[ east ] * velocity_x_east + pressure_east ) )* hxInverse
		+ ( ( rho_u[ north ] * velocity_y_north )