Refactoring the transport equation problem. (aa5b83f8) · Commits · TNL / tnl-dev

examples/advection/LaxFridrichs.h

+66 −31

Original line number	Diff line number	Diff line
		@@ -4,6 +4,7 @@
		#include <TNL/Containers/Vector.h>
		#include <TNL/Meshes/Grid.h>
		#include <TNL/Functions/VectorField.h>
		#include <TNL/SharedPointer.h>

		namespace TNL {

		@@ -26,6 +27,7 @@ class LaxFridrichs< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index,
		public:

		typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType;
		typedef SharedPointer< MeshType > MeshPointer;
		static const int Dimensions = MeshType::getMeshDimensions();
		typedef typename MeshType::CoordinatesType CoordinatesType;
		typedef Real RealType;
		@@ -34,13 +36,25 @@ class LaxFridrichs< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index,
		typedef Functions::MeshFunction< MeshType > MeshFunctionType;
		typedef VelocityFunction VelocityFunctionType;
		typedef Functions::VectorField< Dimensions, VelocityFunctionType > VelocityFieldType;
		typedef SharedPointer< VelocityFieldType > VelocityFieldPointer;

		LaxFridrichs() : artificialViscosity( 1.0 ) {}
		static void configSetup( Config::ConfigDescription& config,
		const String& prefix = "" )
		{
		config.addEntry< double >( "viscosity", "Value of artificial (numerical) viscosity in the Lax-Fridrichs shceme", 1.0 );
		}

		LaxFridrichs()
		: artificialViscosity( 1.0 ) {}

		bool setup( const Config::ParameterContainer& parameters,
		LaxFridrichs( const VelocityFieldPointer& velocityField )
		: artificialViscosity( 1.0 ), velocityField( velocityField ) {}

		bool setup( const MeshPointer& meshPointer,
		const Config::ParameterContainer& parameters,
		const String& prefix = "" )
		{
		return true;
		return this->velocityField->setup( meshPointer, parameters, prefix + "velocity-field-" );
		}

		static String getType();
		@@ -55,12 +69,12 @@ class LaxFridrichs< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index,
		this->tau = tau;
		};

		VelocityFieldType& getVelocityField()
		void setVelocityField( const VelocityFieldPointer& velocityField )
		{
		return this->velocityField;
		this->velocityField = velocityField;
		}

		const VelocityFieldType& getVelocityField() const
		const VelocityFieldPointer& getVelocityField() const
		{
		return this->velocityField;
		}
		@@ -81,7 +95,7 @@ class LaxFridrichs< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index,
		const IndexType& west = neighbourEntities.template getEntityIndex< -1 >();
		typedef Functions::FunctionAdapter< MeshType, VelocityFunctionType > FunctionAdapter;
		return ( 0.5 / this->tau ) * this->artificialViscosity * ( u[ west ]- 2.0 * u[ center ] + u[ east ] ) -
		FunctionAdapter::getValue( this->velocityField[ 0 ], entity, time ) * ( u[ east ] - u[west] ) * hxInverse * 0.5;
		FunctionAdapter::getValue( this->velocityField.getData()[ 0 ], entity, time ) * ( u[ east ] - u[west] ) * hxInverse * 0.5;
		}

		protected:
		@@ -90,7 +104,7 @@ class LaxFridrichs< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index,

		RealType artificialViscosity;

		VelocityFieldType velocityField;
		VelocityFieldPointer velocityField;
		};

		template< typename MeshReal,
		@@ -104,6 +118,7 @@ class LaxFridrichs< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index,
		public:

		typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType;
		typedef SharedPointer< MeshType > MeshPointer;
		static const int Dimensions = MeshType::getMeshDimensions();
		typedef typename MeshType::CoordinatesType CoordinatesType;
		typedef Real RealType;
		@@ -112,14 +127,25 @@ class LaxFridrichs< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index,
		typedef Functions::MeshFunction< MeshType > MeshFunctionType;
		typedef VelocityFunction VelocityFunctionType;
		typedef Functions::VectorField< Dimensions, VelocityFunctionType > VelocityFieldType;
		typedef SharedPointer< VelocityFieldType > VelocityFieldPointer;

		static void configSetup( Config::ConfigDescription& config,
		const String& prefix = "" )
		{
		config.addEntry< double >( "viscosity", "Value of artificial (numerical) viscosity in the Lax-Fridrichs shceme", 1.0 );
		}

		LaxFridrichs()
		: artificialViscosity( 1.0 ) {}

		bool setup( const Config::ParameterContainer& parameters,
		LaxFridrichs( const VelocityFieldPointer& velocityField )
		: artificialViscosity( 1.0 ), velocityField( velocityField ) {}

		bool setup( const MeshPointer& meshPointer,
		const Config::ParameterContainer& parameters,
		const String& prefix = "" )
		{
		return true;
		return this->velocityField->setup( meshPointer, parameters, prefix + "velocity-field-" );
		}

		static String getType();
		@@ -134,12 +160,12 @@ class LaxFridrichs< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index,
		this->tau = tau;
		};

		VelocityFieldType& getVelocityField()
		void setVelocityField( const VelocityFieldPointer& velocityField )
		{
		return this->velocityField;
		this->velocityField = velocityField;
		}

		const VelocityFieldType& getVelocityField() const
		const VelocityFieldPointer& getVelocityField() const
		{
		return this->velocityField;
		}
		@@ -165,8 +191,8 @@ class LaxFridrichs< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index,

		typedef Functions::FunctionAdapter< MeshType, VelocityFunctionType > FunctionAdapter;
		return ( 0.25 / this->tau ) * this->artificialViscosity * ( u[ west ] + u[ east ] + u[ north ] + u[ south ] - 4.0 * u[ center ] ) -
		0.5 * ( FunctionAdapter::getValue( this->velocityField[ 0 ], entity, time ) * ( u[ east ] - u[ west ] ) * hxInverse +
		FunctionAdapter::getValue( this->velocityField[ 1 ], entity, time ) * ( u[ north ] - u[ south ] ) * hyInverse );
		0.5 * ( FunctionAdapter::getValue( this->velocityField.getData()[ 0 ], entity, time ) * ( u[ east ] - u[ west ] ) * hxInverse +
		FunctionAdapter::getValue( this->velocityField.getData()[ 1 ], entity, time ) * ( u[ north ] - u[ south ] ) * hyInverse );
		}

		protected:
		@@ -175,7 +201,7 @@ class LaxFridrichs< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index,

		RealType artificialViscosity;

		VelocityFieldType velocityField;
		VelocityFieldPointer velocityField;
		};

		template< typename MeshReal,
		@@ -189,6 +215,7 @@ class LaxFridrichs< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index,
		public:

		typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;
		typedef SharedPointer< MeshType > MeshPointer;
		static const int Dimensions = MeshType::getMeshDimensions();
		typedef typename MeshType::CoordinatesType CoordinatesType;
		typedef Real RealType;
		@@ -197,14 +224,25 @@ class LaxFridrichs< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index,
		typedef Functions::MeshFunction< MeshType > MeshFunctionType;
		typedef VelocityFunction VelocityFunctionType;
		typedef Functions::VectorField< Dimensions, VelocityFunctionType > VelocityFieldType;
		typedef SharedPointer< VelocityFieldType > VelocityFieldPointer;

		static void configSetup( Config::ConfigDescription& config,
		const String& prefix = "" )
		{
		config.addEntry< double >( "viscosity", "Value of artificial (numerical) viscosity in the Lax-Fridrichs shceme", 1.0 );
		}

		LaxFridrichs()
		: artificialViscosity( 1.0 ) {}

		bool setup( const Config::ParameterContainer& parameters,
		LaxFridrichs( const VelocityFieldPointer& velocityField )
		: artificialViscosity( 1.0 ), velocityField( velocityField ) {}

		bool setup( const MeshPointer& meshPointer,
		const Config::ParameterContainer& parameters,
		const String& prefix = "" )
		{
		return true;
		return this->velocityField->setup( meshPointer, parameters, prefix + "velocity-field-" );
		}

		static String getType();
		@@ -219,12 +257,12 @@ class LaxFridrichs< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index,
		this->tau = tau;
		};

		VelocityFieldType& getVelocityField()
		void setVelocityField( const VelocityFieldPointer& velocityField )
		{
		return this->velocityField;
		this->velocityField = velocityField;
		}

		const VelocityFieldType& getVelocityField() const
		const VelocityFieldPointer& getVelocityField() const
		{
		return this->velocityField;
		}
		@@ -252,9 +290,9 @@ class LaxFridrichs< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index,

		typedef Functions::FunctionAdapter< MeshType, VelocityFunctionType > FunctionAdapter;
		return ( 0.25 / this->tau ) * this->artificialViscosity * ( u[ west ] + u[ east ] + u[ north ] + u[ south ] + u[ up ] + u[ down ] - 6.0 * u[ center ] ) -
		0.5 * ( FunctionAdapter::getValue( this->velocityField[ 0 ], entity, time ) * ( u[ east ] - u[ west ] ) * hxInverse +
		FunctionAdapter::getValue( this->velocityField[ 1 ], entity, time ) * ( u[ north ] - u[ south ] ) * hyInverse +
		FunctionAdapter::getValue( this->velocityField[ 2 ], entity, time ) * ( u[ up ] - u[ down ] ) * hzInverse );
		0.5 * ( FunctionAdapter::getValue( this->velocityField.getData()[ 0 ], entity, time ) * ( u[ east ] - u[ west ] ) * hxInverse +
		FunctionAdapter::getValue( this->velocityField.getData()[ 1 ], entity, time ) * ( u[ north ] - u[ south ] ) * hyInverse +
		FunctionAdapter::getValue( this->velocityField.getData()[ 2 ], entity, time ) * ( u[ up ] - u[ down ] ) * hzInverse );
		}

		protected:
		@@ -263,13 +301,10 @@ class LaxFridrichs< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index,

		RealType artificialViscosity;

		VelocityFieldType velocityField;
		VelocityFieldPointer velocityField;
		};


		} // namespace TNL


		//#include "LaxFridrichs_impl.h"

		#endif /* LaxFridrichs_H */

examples/advection/advection.h

+12 −20

Original line number	Diff line number	Diff line
		@@ -4,6 +4,8 @@
		#include <TNL/Operators/DirichletBoundaryConditions.h>
		#include <TNL/Operators/NeumannBoundaryConditions.h>
		#include <TNL/Functions/Analytic/Constant.h>
		#include <TNL/Functions/VectorField.h>
		#include <TNL/Meshes/Grid.h>
		#include "advectionProblem.h"
		#include "LaxFridrichs.h"
		#include "advectionRhs.h"
		@@ -28,29 +30,19 @@ template< typename ConfigTag >class advectionConfig
		public:
		static void configSetup( Config::ConfigDescription& config )
		{
		config.addDelimiter( "advection settings:" );
		config.addDelimiter( "Advection settings:" );
		config.addEntry< String >( "velocity-field", "Type of velocity field.", "constant" );
		config.addEntryEnum< String >( "constant" );
		//config.addEntryEnum< String >( "file" );
		Functions::VectorField< 3, Functions::Analytic::Constant< 3 > >::configSetup( config, "velocity-field-" );

		typedef Meshes::Grid< 3 > MeshType;
		LaxFridrichs< MeshType >::ConfigSetup( config, "lax-fridrichs" );

		config.addEntry< String >( "boundary-conditions-type", "Choose the boundary conditions type.", "dirichlet");
		config.addEntryEnum< String >( "dirichlet" );
		config.addEntryEnum< String >( "neumann" );
		config.addEntry< double >( "boundary-conditions-constant", "This sets a value in case of the constant boundary conditions." );
		config.addEntry< double >( "artifical-viscosity", "This sets value of artifical viscosity (default 1)", 1.0);
		config.addEntry< String >( "begin", "choose begin type", "sin");
		config.addEntryEnum< String >( "sin");
		config.addEntryEnum< String >( "sin_square");
		config.addEntry< double >( "advection-speedX", "This sets value of advection speed in X direction (default 1)" , 1.0);
		config.addEntry< double >( "advection-speedY", "This sets value of advection speed in Y direction (default 1)" , 1.0);
		config.addEntry< String >( "move", "choose movement type", "advection");
		config.addEntryEnum< String >( "advection");
		config.addEntryEnum< String >( "rotation");
		config.addEntry< int >( "dimension", "choose movement typeproblem dimension", 1);
		config.addEntryEnum< int >( 1 );
		config.addEntryEnum< int >( 2 );
		config.addEntry< double >( "realSize", "Real size of scheme", 1.0);

		/****
		* Add definition of your solver command line arguments.
		*/

		}
		};

examples/advection/advectionProblem_impl.h

+7 −6

Original line number	Diff line number	Diff line
		@@ -55,15 +55,16 @@ setup( const MeshPointer& meshPointer,
		const Config::ParameterContainer& parameters,
		const String& prefix )
		{
		this->velocityType = parameters.getParameter< String >( "move" );
		const double artificalViscosity = parameters.getParameter< double >( "artifical-viscosity" );
		//this->velocityType = parameters.getParameter< String >( "move" );
		/*const double artificalViscosity = parameters.getParameter< double >( "artifical-viscosity" );
		differentialOperatorPointer->setViscosity(artificalViscosity);
		const double advectionSpeedX = parameters.getParameter< double >( "advection-speedX" );
		differentialOperatorPointer->setAdvectionSpeedX(advectionSpeedX);
		const double advectionSpeedY = parameters.getParameter< double >( "advection-speedY" );
		differentialOperatorPointer->setAdvectionSpeedY(advectionSpeedY);
		differentialOperatorPointer->setAdvectionSpeedY(advectionSpeedY);*/

		if( ! this->boundaryConditionPointer->setup( meshPointer, parameters, prefix + "boundary-conditions-" ) \|\|
		if( ! this->differentialOperatorPointer->setup( meshPointer, parameters, prefix + "advection-" ) \|\|
		! this->boundaryConditionPointer->setup( meshPointer, parameters, prefix + "boundary-conditions-" ) \|\|
		! this->rightHandSidePointer->setup( parameters, prefix + "right-hand-side-" ) )
		return false;
		return true;
		@@ -128,7 +129,7 @@ advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >
		setupLinearSystem( const MeshPointer& mesh,
		Matrix& matrix )
		{
		const IndexType dofs = this->getDofs( mesh );
		/*const IndexType dofs = this->getDofs( mesh );
		typedef typename Matrix::ObjectType::CompressedRowsLengthsVector CompressedRowsLengthsVectorType;
		SharedPointer< CompressedRowsLengthsVectorType > rowLengths;
		if( ! rowLengths->setSize( dofs ) )
		@@ -140,7 +141,7 @@ setupLinearSystem( const MeshPointer& mesh,
		rowLengths );
		matrix->setDimensions( dofs, dofs );
		if( ! matrix->setCompressedRowsLengths( *rowLengths ) )
		return false;
		return false;*/
		return true;
		}

src/TNL/Functions/VectorField.h

+42 −7

Original line number	Diff line number	Diff line
		@@ -28,8 +28,26 @@ class VectorField

		typedef Function FunctionType;

		bool setup( const Config::ParameterContainer& parameters,
		const String& prefix = "" );
		static void configSetup( Config::ConfigDescription& config,
		const String& prefix = "" )
		{
		for( int i = 0; i < Dimensions; i++ )
		FunctionType::configSetup( config, prefix + String( i ) + "-" );
		}

		template< typename MeshPointer >
		bool setup( const MeshPointer& meshPointer,
		const Config::ParameterContainer& parameters,
		const String& prefix = "" )
		{
		for( int i = 0; i < Dimensions; i++ )
		if( ! vectorField[ 0 ].setup( meshPointer, parameters, prefix + String( i ) + "-" ) )
		{
		std::cerr << "Unable to setup " << i << "-th coordinate of the vector field." << std::endl;
		return false;
		}
		return true;
		}

		__cuda_callable__
		const FunctionType& operator[]( int i ) const
		@@ -58,10 +76,17 @@ class VectorField< Dimensions, MeshFunction< Mesh, MeshEntityDimensions, Real >
		public:

		typedef Mesh MeshType;
		typedef SharedPointer< MeshType > MeshPointer;
		typedef MeshFunction< MeshType, MeshEntityDimensions, Real > FunctionType;
		typedef typename MeshType::DeviceType DeviceType;
		typedef typename MeshType::IndexType IndexType;
		typedef SharedPointer< MeshType > MeshPointer;

		static void configSetup( Config::ConfigDescription& config,
		const String& prefix = "" )
		{
		for( int i = 0; i < Dimensions; i++ )
		FunctionType::configSetup( config, prefix + String( i ) + "-" );
		}

		VectorField() {};

		@@ -71,8 +96,18 @@ class VectorField< Dimensions, MeshFunction< Mesh, MeshEntityDimensions, Real >
		this->vectorField[ i ].setMesh( meshPointer );
		};

		bool setup( const Config::ParameterContainer& parameters,
		const String& prefix = "" );
		bool setup( const MeshPointer& meshPointer,
		const Config::ParameterContainer& parameters,
		const String& prefix = "" )
		{
		for( int i = 0; i < Dimensions; i++ )
		if( ! vectorField[ 0 ].setup( meshPointer, parameters, prefix + String( i ) + "-" ) )
		{
		std::cerr << "Unable to setup " << i << "-th coordinate of the vector field." << std::endl;
		return false;
		}
		return true;
		}

		__cuda_callable__
		const FunctionType& operator[]( int i ) const