Implementing The Lax Fridrichs scheme for deformed grids. (e42fb244) · Commits · TNL / tnl-dev

examples/navier-stokes/Makefile

+6 −5

Original line number	Diff line number	Diff line
		@@ -9,10 +9,11 @@ INSTALL_DIR = ${HOME}/local
		CXX = g++
		CUDA_CXX = nvcc
		OMP_FLAGS = -DHAVE_OPENMP -fopenmp
		#CXX_FLAGS = -std=gnu++0x -I$(TNL_INCLUDE_DIR) -O0 -g -DDEBUG $(OMP_FLAGS)
		#CXX_FLAGS = -std=gnu++0x -I$(TNL_INCLUDE_DIR) -O3 $(OMP_FLAGS)
		CXX_FLAGS = -std=gnu++0x -I$(TNL_INCLUDE_DIR) -O3 $(OMP_FLAGS) -pg
		#CXX_FLAGS = -DHAVE_NOT_CXX11 -I$(TNL_INCLUDE_DIR) -O3 $(OMP_FLAGS)
		#CXX_FLAGS = -std=gnu++0x -I$(TNL_INCLUDE_DIR) -O0 -g -DDEBUG $(OMP_FLAGS) -DTEMPLATE_EXPLICIT_INSTANTIATION
		CXX_FLAGS = -std=gnu++0x -I$(TNL_INCLUDE_DIR) -O3 $(OMP_FLAGS) -DTEMPLATE_EXPLICIT_INSTANTIATION
		#CXX_FLAGS = -std=gnu++0x -I$(TNL_INCLUDE_DIR) -O3 $(OMP_FLAGS) -pg -DTEMPLATE_EXPLICIT_INSTANTIATION
		#CXX_FLAGS = -DHAVE_NOT_CXX11 -I$(TNL_INCLUDE_DIR) -O3 $(OMP_FLAGS) -DTEMPLATE_EXPLICIT_INSTANTIATION


		LD_FLAGS = -L$(TNL_INSTALL_DIR) -ltnl-0.1 -lgomp

examples/navier-stokes/navierStokesSolver_impl.h

+12 −12

Original line number	Diff line number	Diff line
		@@ -111,7 +111,7 @@ bool navierStokesSolver< Mesh, EulerScheme > :: init( const tnlParameterContaine
		cerr << "Error: y-size must be positive integer number! It is " << meshes. y() << " now." << endl;
		return false;
		}
		this -> mesh. setDimensions( meshes. y(), meshes. x() );
		this -> mesh. setDimensions( meshes. x(), meshes. y() );
		RealType hx = this -> mesh. getParametricStep(). x();
		RealType hy = this -> mesh. getParametricStep(). y();
		mesh. save( tnlString( "mesh.tnl" ) );
		@@ -176,7 +176,7 @@ bool navierStokesSolver< Mesh, EulerScheme > :: setInitialCondition( const tnlPa
		for( IndexType j = 0; j < ySize; j ++ )
		for( IndexType i = 0; i < xSize; i ++ )
		{
		const IndexType c = mesh. getElementIndex( j, i );
		const IndexType c = mesh. getElementIndex( i, j );
		const RealType x = i * hx;
		const RealType y = j * hy;

		@@ -269,7 +269,7 @@ void navierStokesSolver< Mesh, EulerScheme > :: updatePhysicalQuantities( const
		for( IndexType j = 0; j < ySize; j ++ )
		for( IndexType i = 0; i < xSize; i ++ )
		{
		IndexType c = mesh. getElementIndex( j, i );
		IndexType c = mesh. getElementIndex( i, j );
		u1[ c ] = rho_u1[ c ] / rho[ c ];
		u2[ c ] = rho_u2[ c ] / rho[ c ];
		p[ c ] = rho[ c ] * this -> R * this -> T;
		@@ -359,10 +359,10 @@ void navierStokesSolver< Mesh, EulerScheme > :: GetExplicitRHS( const RealType&
		*/
		for( IndexType i = 0; i < xSize; i ++ )
		{
		const IndexType c1 = mesh. getElementIndex( 0, i );
		const IndexType c2 = mesh. getElementIndex( 1, i );
		const IndexType c3 = mesh. getElementIndex( ySize - 1, i );
		const IndexType c4 = mesh. getElementIndex( ySize - 2, i );
		const IndexType c1 = mesh. getElementIndex( i, 0 );
		const IndexType c2 = mesh. getElementIndex( i, 1 );
		const IndexType c3 = mesh. getElementIndex( i, ySize - 1 );
		const IndexType c4 = mesh. getElementIndex( i, ySize - 2 );

		RealType x = i * hx / mesh. getProportions(). x();

		@@ -388,10 +388,10 @@ void navierStokesSolver< Mesh, EulerScheme > :: GetExplicitRHS( const RealType&
		}
		for( IndexType j = 0; j < ySize; j ++ )
		{
		const IndexType c1 = mesh. getElementIndex( j, 0 );
		const IndexType c2 = mesh. getElementIndex( j, 1 );
		const IndexType c3 = mesh. getElementIndex( j, xSize - 1 );
		const IndexType c4 = mesh. getElementIndex( j, xSize - 2 );
		const IndexType c1 = mesh. getElementIndex( 0, j );
		const IndexType c2 = mesh. getElementIndex( 1, j );
		const IndexType c3 = mesh. getElementIndex( xSize - 1, j );
		const IndexType c4 = mesh. getElementIndex( xSize - 2, j );

		RealType y = j * hy / mesh. getProportions(). y();

		@@ -422,7 +422,7 @@ void navierStokesSolver< Mesh, EulerScheme > :: GetExplicitRHS( const RealType&
		for( IndexType j = 0; j < ySize; j ++ )
		for( IndexType i = 0; i < xSize; i ++ )
		{
		IndexType c = mesh. getElementIndex( j, i );
		IndexType c = mesh. getElementIndex( i, j );
		if( i == 0 \|\| j == 0 \|\|
		i == xSize - 1 \|\| j == ySize - 1 )
		{

src/implementation/mesh/tnlGrid2D_impl.h

+35 −11

Original line number	Diff line number	Diff line
		@@ -58,7 +58,7 @@ template< typename Real,
		typename Device,
		typename Index,
		template< int, typename, typename, typename > class Geometry >
		void tnlGrid< 2, Real, Device, Index, Geometry > :: setDimensions( const Index ySize, const Index xSize )
		void tnlGrid< 2, Real, Device, Index, Geometry > :: setDimensions( const Index xSize, const Index ySize )
		{
		tnlAssert( xSize > 0,
		cerr << "The number of Elements along x-axis must be larger than 0." );
		@@ -80,7 +80,7 @@ template< typename Real,
		template< int, typename, typename, typename > class Geometry >
		void tnlGrid< 2, Real, Device, Index, Geometry > :: setDimensions( const tnlTuple< 2, Index >& dimensions )
		{
		this -> setDimensions( dimensions. y(), dimensions. x() );
		this -> setDimensions( dimensions. x(), dimensions. y() );
		}

		template< typename Real,
		@@ -169,7 +169,7 @@ template< typename Real,
		typename Device,
		typename Index,
		template< int, typename, typename, typename > class Geometry >
		Index tnlGrid< 2, Real, Device, Index, Geometry > :: getElementIndex( const Index j, const Index i ) const
		Index tnlGrid< 2, Real, Device, Index, Geometry > :: getElementIndex( const Index i, const Index j ) const
		{
		tnlAssert( i < dimensions. x(),
		cerr << "Index i ( " << i
		@@ -183,6 +183,21 @@ Index tnlGrid< 2, Real, Device, Index, Geometry > :: getElementIndex( const Inde
		return j * this -> dimensions. x() + i;
		}

		template< typename Real,
		typename Device,
		typename Index,
		template< int, typename, typename, typename > class Geometry >
		void tnlGrid< 2, Real, Device, Index, Geometry > :: getElementCoordinates( const Index element,
		tnlTuple< 2, Index >& coordinates ) const
		{
		tnlAssert( i >= 0 && i < dofs,
		cerr << " i = " << i << " dofs = " << dofs
		<< " in tnlGrid " << this -> getName(); );

		coordinates. x() = element % this -> dimensions. x();
		coordinates. y() = element / this -> dimensions. x();
		}

		template< typename Real,
		typename Device,
		typename Index,
		@@ -213,16 +228,25 @@ template< typename Real,
		typename Device,
		typename Index,
		template< int, typename, typename, typename > class Geometry >
		Real tnlGrid< 2, Real, Device, Index, Geometry > :: getElementMeasure( const Index j,
		const Index i ) const
		void tnlGrid< 2, Real, Device, Index, Geometry > :: getElementCenter( const tnlTuple< 2, Index >& coordinates,
		tnlTuple< 2, Real >& center ) const
		{
		return geometry. getElementMeasure( j, i );
		geometry. getElementCenter( origin, coordinates, center );
		}

		template< typename Real,
		typename Device,
		typename Index,
		template< int, typename, typename, typename > class Geometry >
		Real tnlGrid< 2, Real, Device, Index, Geometry > :: getElementMeasure( const tnlTuple< 2, Index >& coordinates ) const
		{
		return geometry. getElementMeasure( coordinates );
		}

		/*template< typename Real,
		typename Device,
		typename Index,
		template< int, typename, typename, typename > class Geometry >
		template< Index dy, Index dx >
		Real tnlGrid< 2, Real, Device, Index, Geometry > :: getElementsDistance( const Index j,
		const Index i ) const
		@@ -250,7 +274,7 @@ tnlTuple< 2, Real > tnlGrid< 2, Real, Device, Index, Geometry > :: getEdgeNormal
		const Index i ) const
		{
		return geometry. getEdgeNormal< dy, dx >( j, i );
		}
		}*/

		template< typename Real,
		typename Device,
		@@ -349,7 +373,7 @@ bool tnlGrid< 2, Real, Device, Index, Geometry > :: write( const MeshFunction& f
		{
		const RealType x = this -> getOrigin(). x() + i * hx;
		const RealType y = this -> getOrigin(). y() + j * hy;
		file << x << " " << " " << y << " " << function[ this -> getElementIndex( j, i ) ] << endl;
		file << x << " " << " " << y << " " << function[ this -> getElementIndex( i, j ) ] << endl;
		}
		file << endl;
		}

src/implementation/mesh/tnlGrid3D_impl.h

+4 −12

Original line number	Diff line number	Diff line
		@@ -55,7 +55,7 @@ template< typename Real,
		typename Device,
		typename Index,
		template< int, typename, typename, typename > class Geometry >
		void tnlGrid< 3, Real, Device, Index, Geometry > :: setDimensions( const Index zSize, const Index ySize, const Index xSize )
		void tnlGrid< 3, Real, Device, Index, Geometry > :: setDimensions( const Index xSize, const Index ySize, const Index zSize )
		{
		tnlAssert( xSize > 1,
		cerr << "The number of Elements along x-axis must be larger than 1." );
		@@ -76,15 +76,7 @@ template< typename Real,
		template< int, typename, typename, typename > class Geometry >
		void tnlGrid< 3, Real, Device, Index, Geometry > :: setDimensions( const tnlTuple< 3, Index >& dimensions )
		{
		tnlAssert( dimensions. x() > 1,
		cerr << "The number of Elements along x-axis must be larger than 1." );
		tnlAssert( dimensions. y() > 1,
		cerr << "The number of Elements along y-axis must be larger than 1." );
		tnlAssert( dimensions. z() > 1,
		cerr << "The number of Elements along z-axis must be larger than 1." );

		this -> dimensions = dimensions;
		dofs = this -> dimensions. x() * this -> dimensions. y() * this -> dimensions. z();
		this -> setDimensions( this -> dimensions. x(), this -> dimensions. y(), this -> dimensions. z() );
		}

		template< typename Real,
		@@ -177,7 +169,7 @@ template< typename Real,
		typename Device,
		typename Index,
		template< int, typename, typename, typename > class Geometry >
		Index tnlGrid< 3, Real, Device, Index, Geometry > :: getElementIndex( const Index k, const Index j, const Index i ) const
		Index tnlGrid< 3, Real, Device, Index, Geometry > :: getElementIndex( const Index i, const Index j, const Index k ) const
		{
		tnlAssert( i < dimensions. x(),
		cerr << "Index i ( " << i
		@@ -297,7 +289,7 @@ bool tnlGrid< 3, Real, Device, Index, Geometry > :: write( const MeshFunction& f
		{
		const RealType x = this -> getOrigin(). x() + i * hx;
		const RealType y = this -> getOrigin(). y() + j * hy;
		//file << x << " " << " " << y << " " << function[ this -> getElementIndex( j, i ) ] << endl;
		//file << x << " " << " " << y << " " << function[ this -> getElementIndex( i, j ) ] << endl;
		}
		file << endl;
		}

src/implementation/mesh/tnlIdenticalGridGeometry_impl.h

+22 −24

Original line number	Diff line number	Diff line
		@@ -45,11 +45,11 @@ const tnlTuple< 1, Real >& tnlIdenticalGridGeometry< 1, Real, Device, Index > ::
		template< typename Real,
		typename Device,
		typename Index >
		void tnlIdenticalGridGeometry< 1, Real, Device, Index > :: getElementCoordinates( const Index i,
		const tnlTuple< 1, Real >& origin,
		tnlTuple< 1, Real >& coordinates ) const
		void tnlIdenticalGridGeometry< 1, Real, Device, Index > :: getElementCenter( const tnlTuple< 1, Real >& origin,
		const tnlTuple< 1, Index >& coordinates,
		tnlTuple< 1, Real >& center ) const
		{
		coordinates. x() = ( i + 0.5 ) * parametricStep. x();
		center. x() = ( coordinates. x() + 0.5 ) * parametricStep. x();
		}

		template< typename Real,
		@@ -145,20 +145,18 @@ const tnlTuple< 2, Real >& tnlIdenticalGridGeometry< 2, Real, Device, Index > ::
		template< typename Real,
		typename Device,
		typename Index >
		void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementCoordinates( const Index j,
		const Index i,
		const tnlTuple< 2, Real >& origin,
		tnlTuple< 2, Real >& coordinates ) const
		void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementCenter( const tnlTuple< 2, Real >& origin,
		const tnlTuple< 2, Index >& coordinates,
		tnlTuple< 2, Real >& center ) const
		{
		coordinates. x() = ( i + 0.5 ) * parametricStep. x();
		coordinates. y() = ( j + 0.5 ) * parametricStep. y();
		center. x() = ( coordinates. x() + 0.5 ) * parametricStep. x();
		center. y() = ( coordinates. y() + 0.5 ) * parametricStep. y();
		}

		template< typename Real,
		typename Device,
		typename Index >
		Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementMeasure( const Index j,
		const Index i ) const
		Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementMeasure( const tnlTuple< 2, Index >& coordinates ) const
		{
		return elementMeasure;
		}
		@@ -166,9 +164,9 @@ Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementMeasure( co
		template< typename Real,
		typename Device,
		typename Index >
		template< Index dy, Index dx >
		Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementsDistance( const Index j,
		const Index i ) const
		template< Index dx, Index dy >
		Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementsDistance( const Index i,
		const Index j ) const
		{
		if( dy == 0 && dx == 1 )
		return parametricStep. x();
		@@ -182,9 +180,9 @@ Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getElementsDistance(
		template< typename Real,
		typename Device,
		typename Index >
		template< Index dy, Index dx >
		void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getEdgeCoordinates( const Index j,
		const Index i,
		template< Index dx, Index dy >
		void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getEdgeCoordinates( const Index i,
		const Index j,
		const tnlTuple< 2, Real >& origin,
		tnlTuple< 2, Real >& coordinates ) const
		{
		@@ -195,9 +193,9 @@ void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getEdgeCoordinates( c
		template< typename Real,
		typename Device,
		typename Index >
		template< Index dy, Index dx >
		Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getEdgeLength( const Index j,
		const Index i ) const
		template< Index dx, Index dy >
		Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getEdgeLength( const Index i,
		const Index j ) const
		{
		if( dy == 0 && dx == 1 )
		return parametricStep. y();
		@@ -209,9 +207,9 @@ Real tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getEdgeLength( const
		template< typename Real,
		typename Device,
		typename Index >
		template< Index dy, Index dx >
		void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getEdgeNormal( const Index j,
		const Index i,
		template< Index dx, Index dy >
		void tnlIdenticalGridGeometry< 2, Real, Device, Index > :: getEdgeNormal( const Index i,
		const Index j,
		tnlTuple< 2, Real >& normal ) const
		{
		tnlAssert( ( dx == 0 \|\| dx == 1 \|\| dx == -1 \|\|