Refactoring the Navier-Stokes solver example. (697ee489) · Commits · TNL / tnl-dev

examples/navier-stokes/navierStokesSolver.h

+6 −1

Original line number	Diff line number	Diff line
		@@ -29,6 +29,7 @@
		#include <schemes/gradient/tnlCentralFDMGradient.h>
		#include <schemes/diffusion/tnlLinearDiffusion.h>
		#include <mesh/tnlLinearGridGeometry.h>
		#include <schemes/navier-stokes/tnlNavierStokes.h>

		#include "navierStokesSolverMonitor.h"
		#include "navierStokesBoundaryConditions.h"
		@@ -97,7 +98,7 @@ class navierStokesSolver

		MeshType mesh;

		tnlVector< RealType, DeviceType, IndexType > rho, u1, u2, p;
		//tnlVector< RealType, DeviceType, IndexType > rho, u1, u2, p;

		DofVectorType dofVector, rhsDofVector;

		@@ -105,6 +106,10 @@ class navierStokesSolver

		EulerScheme eulerScheme;

		tnlNavierStokes< EulerScheme,
		tnlLinearDiffusion< MeshType >,
		navierStokesBoundaryConditions< MeshType > > navierStokesScheme;

		tnlLinearDiffusion< MeshType > u1Viscosity, u2Viscosity;

		tnlCentralFDMGradient< MeshType > pressureGradient;

examples/navier-stokes/navierStokesSolver_impl.h

+25 −26

Original line number	Diff line number	Diff line
		@@ -53,10 +53,7 @@ navierStokesSolver< Mesh, EulerScheme > :: navierStokesSolver()
		p_0( 0.0 ),
		gravity( 0.0 )
		{
		this -> rho. setName( "rho-g" );
		this -> u1. setName( "u1-g");
		this -> u2. setName( "u2-g" );
		this -> p. setName( "p" );

		this -> mesh. setName( "navier-stokes-mesh" );
		this -> dofVector. setName( "navier-stokes-dof-vector" );
		}
		@@ -151,11 +148,6 @@ bool navierStokesSolver< Mesh, EulerScheme > :: init( const tnlParameterContaine
		*/
		const IndexType variablesNumber = 3;

		rho.setSize( mesh.getDofs() );
		u1.setSize( mesh.getDofs() );
		u2.setSize( mesh.getDofs() );
		p.setSize( mesh.getDofs() );

		dofVector. setSize( variablesNumber * mesh. getDofs() );
		rhsDofVector. setLike( dofVector );
		this->boundaryConditions.setMesh( this->mesh );
		@@ -163,14 +155,17 @@ bool navierStokesSolver< Mesh, EulerScheme > :: init( const tnlParameterContaine
		/****
		* Set-up numerical scheme
		*/
		pressureGradient. setFunction( p );
		pressureGradient.setFunction( navierStokesScheme.getPressure() );
		pressureGradient.bindMesh( this -> mesh );
		this->eulerScheme. bindMesh( this -> mesh );
		this->eulerScheme. setPressureGradient( this -> pressureGradient );
		this->u1Viscosity. bindMesh( this -> mesh );
		this -> u1Viscosity. setFunction( this -> u1 );
		this->u1Viscosity. setFunction( this -> navierStokesScheme.getU1() );
		this->u2Viscosity. bindMesh( this -> mesh );
		this -> u2Viscosity. setFunction( this -> u2 );
		this->u2Viscosity. setFunction( this -> navierStokesScheme.getU2() );
		navierStokesScheme.setAdvectionScheme( this->eulerScheme );
		navierStokesScheme.setMesh( this->mesh );
		//navierStokesScheme.setDifusionScheme( this)
		return true;
		}

		@@ -229,11 +224,15 @@ bool navierStokesSolver< Mesh, EulerScheme > :: makeSnapshot( const RealType& t,
		rho_u2_t. bind( & this->rhsDofVector.getData()[ 2 * dofs ], dofs );


		updatePhysicalQuantities( rho, rho_u1, rho_u2 );
		navierStokesScheme.updatePhysicalQuantities( rho, rho_u1, rho_u2 );
		tnlVector< RealType, DeviceType, IndexType > u;
		u. setLike( u1 );
		u. setLike( navierStokesScheme.getU1() );
		for( IndexType i = 0; i < this -> u1. getSize(); i ++ )
		u[ i ] = sqrt( u1[ i ] * u1[ i ] + u2[ i ] * u2[ i ] );
		{
		const RealType& u1 = navierStokesScheme.getU1()[ i ];
		const RealType& u2 = navierStokesScheme.getU2()[ i ];
		u[ i ] = sqrt( u1 * u1 + u2 * u2 );
		}
		tnlString fileName;
		/*FileNameBaseNumberEnding( "u-1-", step, 5, ".tnl", fileName );
		if( ! this -> u1. save( fileName ) )
		@@ -247,7 +246,7 @@ bool navierStokesSolver< Mesh, EulerScheme > :: makeSnapshot( const RealType& t,
		FileNameBaseNumberEnding( "u-", step, 5, ".tnl", fileName );
		if( ! u. save( fileName ) )
		return false;
		FileNameBaseNumberEnding( "rho-t-", step, 5, ".tnl", fileName );
		/*FileNameBaseNumberEnding( "rho-t-", step, 5, ".tnl", fileName );
		if( ! rho_t. save( fileName ) )
		return false;
		FileNameBaseNumberEnding( "rho-u1-t-", step, 5, ".tnl", fileName );
		@@ -256,10 +255,10 @@ bool navierStokesSolver< Mesh, EulerScheme > :: makeSnapshot( const RealType& t,

		FileNameBaseNumberEnding( "rho-u2-t-", step, 5, ".tnl", fileName );
		if( ! rho_u2_t. save( fileName ) )
		return false;
		return false;*/

		FileNameBaseNumberEnding( "rho-", step, 5, ".tnl", fileName );
		if( ! rho. save( fileName ) )
		if( ! navierStokesScheme.getRho(). save( fileName ) )
		return false;
		FileNameBaseNumberEnding( "rho-u1-", step, 5, ".tnl", fileName );
		if( ! rho_u1. save( fileName ) )
		@@ -271,7 +270,7 @@ bool navierStokesSolver< Mesh, EulerScheme > :: makeSnapshot( const RealType& t,
		return true;
		}

		template< typename Mesh, typename EulerScheme >
		/*template< typename Mesh, typename EulerScheme >
		template< typename Vector >
		void navierStokesSolver< Mesh, EulerScheme > :: updatePhysicalQuantities( const Vector& dofs_rho,
		const Vector& rho_u1,
		@@ -295,7 +294,7 @@ void navierStokesSolver< Mesh, EulerScheme > :: updatePhysicalQuantities( const
		this->p[ c ] = this->rho[ c ] * this -> R * this -> T;
		}
		}
		}
		}*/

		template< typename Mesh, typename EulerScheme >
		void navierStokesSolver< Mesh, EulerScheme > :: GetExplicitRHS( const RealType& time,
		@@ -319,7 +318,7 @@ void navierStokesSolver< Mesh, EulerScheme > :: GetExplicitRHS( const RealType&
		rho_u1_t. bind( & fu. getData()[ dofs ], dofs );
		rho_u2_t. bind( & fu. getData()[ 2 * dofs ], dofs );

		updatePhysicalQuantities( dofs_rho, rho_u1, rho_u2 );
		navierStokesScheme.updatePhysicalQuantities( dofs_rho, rho_u1, rho_u2 );

		/****
		* Compute RHS

src/implementation/schemes/navier-stokes/tnlNavierStokes_impl.h

+152 −44

Original line number	Diff line number	Diff line
		@@ -21,76 +21,181 @@
		#include <schemes/navier-stokes/tnlNavierStokes.h>

		template< typename AdvectionScheme,
		typename DiffusionScheme >
		tnlNavierStokes< AdvectionScheme, DiffusionScheme >::tnlNavierStokes()
		typename DiffusionScheme,
		typename BoundaryConditions >
		tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::tnlNavierStokes()
		: advection( 0 ),
		diffusion( 0 )
		{
		this->rho.setName( "navier-stokes-rho" );
		this->u1.setName( "navier-stokes-u1");
		this->u2.setName( "navier-stokes-u2" );
		this->p.setName( "navier-stokes-p" );
		}

		template< typename AdvectionScheme,
		typename DiffusionScheme >
		tnlString tnlNavierStokes< AdvectionScheme, DiffusionScheme >::getTypeStatic()
		typename DiffusionScheme,
		typename BoundaryConditions >
		tnlString tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::getTypeStatic()
		{
		return tnlString( "tnlNavierStokes< " ) +
		AdvectionScheme::getTypeStatic() + ", "
		AdvectionScheme::getTypeStatic() + ", " +
		DiffusionScheme::getTypeStatic() + " >";
		}

		template< typename AdvectionScheme,
		typename DiffusionScheme >
		void tnlNavierStokes< AdvectionScheme, DiffusionScheme >::setAdvectionScheme( AdvectionSchemeType& advection )
		typename DiffusionScheme,
		typename BoundaryConditions >
		void tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::setAdvectionScheme( AdvectionSchemeType& advection )
		{
		this->advection = advection;
		this->advection = &advection;
		}

		template< typename AdvectionScheme,
		typename DiffusionScheme >
		void tnlNavierStokes< AdvectionScheme, DiffusionScheme >::setDiffusionScheme( DiffusionSchemeType& diffusion )
		typename DiffusionScheme,
		typename BoundaryConditions >
		void tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::setDiffusionScheme( DiffusionSchemeType& diffusion )
		{
		this->diffusion = diffusion;
		this->diffusion = &diffusion;
		}

		template< typename AdvectionScheme,
		typename DiffusionScheme >
		void tnlNavierStokes< AdvectionScheme, DiffusionScheme >::setBoundaryConditions( BoundaryConditionsType& boundaryConditions )
		typename DiffusionScheme,
		typename BoundaryConditions >
		void tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::setBoundaryConditions( BoundaryConditionsType& boundaryConditions )
		{
		this->boundaryConditions = boundaryConditions;
		this->boundaryConditions = &boundaryConditions;
		}

		template< typename AdvectionScheme,
		typename DiffusionScheme >
		void tnlNavierStokes< AdvectionScheme, DiffusionScheme > :: updatePhysicalQuantities( const Vector& rho,
		const Vector& rho_u1,
		const Vector& rho_u2 )
		typename DiffusionScheme,
		typename BoundaryConditions >
		void tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::setMesh( MeshType& mesh )
		{
		this->mesh = &mesh;
		this->rho.setSize( this->mesh->getDofs() );
		this->u1.setSize( this->mesh->getDofs() );
		this->u2.setSize( this->mesh->getDofs() );
		this->p.setSize( this->mesh->getDofs() );
		}

		template< typename AdvectionScheme,
		typename DiffusionScheme,
		typename BoundaryConditions >
		typename tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::VectorType&
		tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::getRho()
		{
		return this->rho;
		}

		template< typename AdvectionScheme,
		typename DiffusionScheme,
		typename BoundaryConditions >
		const typename tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::VectorType&
		tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::getRho() const
		{
		return this->rho;
		}

		template< typename AdvectionScheme,
		typename DiffusionScheme,
		typename BoundaryConditions >
		typename tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::VectorType&
		tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::getU1()
		{
		return this->u1;
		}

		template< typename AdvectionScheme,
		typename DiffusionScheme,
		typename BoundaryConditions >
		const typename tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::VectorType&
		tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::getU1() const
		{
		return this->u1;
		}

		template< typename AdvectionScheme,
		typename DiffusionScheme,
		typename BoundaryConditions >
		typename tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::VectorType&
		tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::getU2()
		{
		return this->u2;
		}

		template< typename AdvectionScheme,
		typename DiffusionScheme,
		typename BoundaryConditions >
		const typename tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::VectorType&
		tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::getU2() const
		{
		return this->u2;
		}

		template< typename AdvectionScheme,
		typename DiffusionScheme,
		typename BoundaryConditions >
		typename tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::VectorType&
		tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::getPressure()
		{
		return this->p;
		}

		template< typename AdvectionScheme,
		typename DiffusionScheme,
		typename BoundaryConditions >
		const typename tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::VectorType&
		tnlNavierStokes< AdvectionScheme, DiffusionScheme, BoundaryConditions >::getPressure() const
		{
		return this->p;
		}

		template< typename AdvectionScheme,
		typename DiffusionScheme,
		typename BoundaryConditions >
		template< typename Vector >
		void tnlNavierStokes< AdvectionScheme,
		DiffusionScheme,
		BoundaryConditions > :: updatePhysicalQuantities( const Vector& dofs_rho,
		const Vector& dofs_rho_u1,
		const Vector& dofs_rho_u2 )
		{
		if( DeviceType :: getDevice() == tnlHostDevice )
		{
		const IndexType& xSize = mesh. getDimensions(). x();
		const IndexType& ySize = mesh. getDimensions(). y();
		//const IndexType& xSize = mesh. getDimensions(). x();
		//const IndexType& ySize = mesh. getDimensions(). y();
		const IndexType size = dofs_rho.getSize();

		#ifdef HAVE_OPENMP
		#pragma omp parallel for
		#endif
		for( IndexType j = 0; j < ySize; j ++ )
		for( IndexType i = 0; i < xSize; i ++ )
		//for( IndexType j = 0; j < ySize; j ++ )
		// for( IndexType i = 0; i < xSize; i ++ )
		for( IndexType c = 0; c < size; c++ )
		{
		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;
		//IndexType c = mesh. getElementIndex( i, j );
		this->rho[ c ] = dofs_rho[ c ];
		this->u1[ c ] = dofs_rho_u1[ c ] / dofs_rho[ c ];
		this->u2[ c ] = dofs_rho_u2[ c ] / dofs_rho[ c ];
		this->p[ c ] = dofs_rho[ c ] * this -> R * this -> T;
		}
		}
		}

		void tnlNavierStokes< AdvectionScheme, DiffusionScheme >::getExplicitRhs( const RealType& time,
		template< typename AdvectionScheme,
		typename DiffusionScheme,
		typename BoundaryConditions >
		void tnlNavierStokes< AdvectionScheme,
		DiffusionScheme,
		BoundaryConditions >::getExplicitRhs( const RealType& time,
		const RealType& tau,
		DofVectorType& u,
		DofVectorType& fu )
		DofVectorType& fu ) const
		{
		tnlAssert( this->advection );
		tnlAssert( this->diffusion );
		tnlAssert( this->boundaryConditions );
		tnlAssert( this->advection, );
		tnlAssert( this->diffusion, );
		tnlAssert( this->boundaryConditions, );

		tnlSharedVector< RealType, DeviceType, IndexType > rho, rho_u1, rho_u2,
		rho_t, rho_u1_t, rho_u2_t;
		@@ -100,9 +205,9 @@ void tnlNavierStokes< AdvectionScheme, DiffusionScheme >::getExplicitRhs( const
		rho_u1. bind( & u. getData()[ dofs ], dofs );
		rho_u2. bind( & u. getData()[ 2 * dofs ], dofs );

		eulerScheme. setRho( rho );
		eulerScheme. setRhoU1( rho_u1 );
		eulerScheme. setRhoU2( rho_u2 );
		advection.setRho( rho );
		advection.setRhoU1( rho_u1 );
		advection.setRhoU2( rho_u2 );

		rho_t.bind( & fu. getData()[ 0 ], dofs );
		rho_u1_t.bind( & fu. getData()[ dofs ], dofs );
		@@ -110,6 +215,9 @@ void tnlNavierStokes< AdvectionScheme, DiffusionScheme >::getExplicitRhs( const

		updatePhysicalQuantities( rho, rho_u1, rho_u2 );

		const IndexType& xSize = this->mesh->getDimensions().x();
		const IndexType& ySize = this->mesh->getDimensions().y();

		#ifdef HAVE_OPENMP
		#pragma omp parallel for
		#endif
		@@ -124,7 +232,7 @@ void tnlNavierStokes< AdvectionScheme, DiffusionScheme >::getExplicitRhs( const
		continue;
		}

		eulerScheme. getExplicitRhs( c,
		this->advection->getExplicitRhs( c,
		rho_t[ c ],
		rho_u1_t[ c ],
		rho_u2_t[ c ],

src/schemes/navier-stokes/tnlNavierStokes.h

+30 −3

Original line number	Diff line number	Diff line
		@@ -19,6 +19,7 @@
		#define TNLNAVIERSTOKES_H_

		#include <core/tnlString.h>
		#include <core/vectors/tnlVector.h>

		template< typename AdvectionScheme,
		typename DiffusionScheme,
		@@ -31,9 +32,11 @@ class tnlNavierStokes
		typedef DiffusionScheme DiffusionSchemeType;
		typedef BoundaryConditions BoundaryConditionsType;
		typedef typename AdvectionScheme::MeshType MeshType;
		typedef typename AdvectionScheme::Real RealType;
		typedef typename AdvectionScheme::Device DeviceType;
		typedef typename AdvectionScheme::Index IndexType;
		typedef typename AdvectionScheme::RealType RealType;
		typedef typename AdvectionScheme::DeviceType DeviceType;
		typedef typename AdvectionScheme::IndexType IndexType;
		typedef tnlVector< RealType, DeviceType, IndexType > VectorType;
		typedef VectorType DofVectorType;

		tnlNavierStokes();

		@@ -45,6 +48,25 @@ class tnlNavierStokes

		void setBoundaryConditions( BoundaryConditionsType& boundaryConditions );

		void setMesh( MeshType& mesh );

		VectorType& getRho();

		const VectorType& getRho() const;

		VectorType& getU1();

		const VectorType& getU1() const;

		VectorType& getU2();

		const VectorType& getU2() const;

		VectorType& getPressure();

		const VectorType& getPressure() const;

		template< typename Vector >
		void updatePhysicalQuantities( const Vector& rho,
		const Vector& rho_u1,
		const Vector& rho_u2 );
		@@ -63,7 +85,12 @@ class tnlNavierStokes

		BoundaryConditionsType* boundaryConditions;

		MeshType* mesh;

		VectorType rho, u1, u2, p;

		};

		#include <implementation/schemes/navier-stokes/tnlNavierStokes_impl.h>

		#endif /* TNLNAVIERSTOKES_H_ */