Adding energy equation to the Navier-Stokes solver. (43c75758) · Commits · TNL / tnl-dev

examples/navier-stokes/navierStokesSolver.h

+2 −7

Original line number	Diff line number	Diff line
		@@ -87,12 +87,7 @@ class navierStokesSolver

		protected:

		RealType regularize( const RealType& r ) const;

		template< typename Vector >
		void updatePhysicalQuantities( const Vector& rho,
		const Vector& rho_u1,
		const Vector& rho_u2 );
		//RealType regularize( const RealType& r ) const;

		ProblemType problem;

		@@ -108,7 +103,7 @@ class navierStokesSolver

		tnlNavierStokesSolver< EulerScheme,
		tnlLinearDiffusion< MeshType >,
		navierStokesBoundaryConditions< MeshType > > navierStokesScheme;
		navierStokesBoundaryConditions< MeshType > > nsSolver;

		tnlLinearDiffusion< MeshType > u1Viscosity, u2Viscosity;

examples/navier-stokes/navierStokesSolver_impl.h

+25 −24

Original line number	Diff line number	Diff line
		@@ -127,46 +127,47 @@ bool navierStokesSolver< Mesh, EulerScheme > :: init( const tnlParameterContaine
		RealType hy = this -> mesh. getParametricStep(). y();
		mesh. refresh();
		mesh. save( tnlString( "mesh.tnl" ) );
		nsSolver.setMesh( this->mesh );

		/****
		* Set-up model coefficients
		*/
		this->p_0 = parameters. GetParameter< double >( "p0" );
		navierStokesScheme.setMu( parameters. GetParameter< double >( "mu") );
		navierStokesScheme.setT( parameters. GetParameter< double >( "T") );
		navierStokesScheme.setR( parameters. GetParameter< double >( "R") );
		navierStokesScheme.setGravity( parameters. GetParameter< double >( "gravity") );
		nsSolver.setMu( parameters. GetParameter< double >( "mu") );
		nsSolver.setT( parameters. GetParameter< double >( "T") );
		nsSolver.setR( parameters. GetParameter< double >( "R") );
		nsSolver.setGravity( parameters. GetParameter< double >( "gravity") );
		if( ! this->boundaryConditions.init( parameters ) )
		return false;

		/****
		* Set-up grid functions
		*/
		const IndexType variablesNumber = 3;

		dofVector. setSize( variablesNumber * mesh. getDofs() );
		dofVector. setSize( nsSolver. getDofs() );
		rhsDofVector. setLike( dofVector );
		nsSolver.bindDofVector( dofVector.getData() );

		/****
		* Set-up boundary conditions
		*/
		this->boundaryConditions.setMesh( this->mesh );
		nsSolver.setBoundaryConditions( this->boundaryConditions );

		/****
		* Set-up numerical scheme
		*/
		navierStokesScheme.setMesh( this->mesh );
		pressureGradient.setFunction( navierStokesScheme.getPressure() );

		pressureGradient.setFunction( nsSolver.getPressure() );
		pressureGradient.bindMesh( this -> mesh );
		this->eulerScheme. bindMesh( this -> mesh );
		this->eulerScheme. setPressureGradient( this -> pressureGradient );
		this->u1Viscosity. bindMesh( this -> mesh );
		this->u1Viscosity. setFunction( this -> navierStokesScheme.getU1() );
		this->u1Viscosity. setFunction( this -> nsSolver.getU1() );
		this->u2Viscosity. bindMesh( this -> mesh );
		this->u2Viscosity. setFunction( this -> navierStokesScheme.getU2() );
		navierStokesScheme.setAdvectionScheme( this->eulerScheme );
		navierStokesScheme.setDiffusionScheme( this->u1Viscosity,
		this->u2Viscosity. setFunction( this -> nsSolver.getU2() );
		nsSolver.setAdvectionScheme( this->eulerScheme );
		nsSolver.setDiffusionScheme( this->u1Viscosity,
		this->u2Viscosity );
		navierStokesScheme.setBoundaryConditions( this->boundaryConditions );
		navierStokesScheme.bindDofVector( dofVector.getData() );

		//navierStokesScheme.setDifusionScheme( this)
		return true;
		}

		@@ -179,8 +180,8 @@ bool navierStokesSolver< Mesh, EulerScheme > :: setInitialCondition( const tnlPa
		rho_u1. bind( & dofVector. getData()[ dofs ], dofs );
		rho_u2. bind( & dofVector. getData()[ 2 * dofs ], dofs );

		dofs_rho. setValue( p_0 / ( this->navierStokesScheme.getR() *
		this->navierStokesScheme.getT() ) );
		dofs_rho. setValue( p_0 / ( this->nsSolver.getR() *
		this->nsSolver.getT() ) );
		rho_u1. setValue( 0.0 );
		rho_u2. setValue( 0.0 );

		@@ -196,8 +197,8 @@ bool navierStokesSolver< Mesh, EulerScheme > :: setInitialCondition( const tnlPa
		const RealType x = i * hx;
		const RealType y = j * hy;

		dofs_rho. setElement( c, p_0 / ( this->navierStokesScheme.getR() *
		this->navierStokesScheme.getT() ) );
		dofs_rho. setElement( c, p_0 / ( this->nsSolver.getR() *
		this->nsSolver.getT() ) );
		rho_u1. setElement( c, 0.0 );
		rho_u2. setElement( c, 0.0 );

		@@ -216,9 +217,9 @@ bool navierStokesSolver< Mesh, EulerScheme > :: makeSnapshot( const RealType& t,
		const IndexType step )
		{
		cout << endl << "Writing output at time " << t << " step " << step << "." << endl;
		if( !navierStokesScheme.writePhysicalVariables( t, step ) )
		if( !nsSolver.writePhysicalVariables( t, step ) )
		return false;
		if( !navierStokesScheme.writeConservativeVariables( t, step ) )
		if( !nsSolver.writeConservativeVariables( t, step ) )
		return false;
		return true;
		}
		@@ -229,7 +230,7 @@ void navierStokesSolver< Mesh, EulerScheme > :: GetExplicitRHS( const RealType&
		DofVectorType& u,
		DofVectorType& fu )
		{
		navierStokesScheme.getExplicitRhs( time, tau, u, fu );
		nsSolver.getExplicitRhs( time, tau, u, fu );
		}

		template< typename Mesh, typename EulerScheme >

src/implementation/solvers/cfd/navier-stokes/tnlNavierStokesSolver_impl.h

+55 −19

Original line number	Diff line number	Diff line
		@@ -36,6 +36,7 @@ tnlNavierStokesSolver< AdvectionScheme, DiffusionScheme, BoundaryConditions >::t
		this->u1.setName( "navier-stokes-u1");
		this->u2.setName( "navier-stokes-u2" );
		this->p.setName( "navier-stokes-p" );
		this->e.setName( "navier-stokes-e" );
		}

		template< typename AdvectionScheme,
		@@ -84,6 +85,7 @@ void tnlNavierStokesSolver< AdvectionScheme, DiffusionScheme, BoundaryConditions
		this->u1.setSize( this->mesh->getDofs() );
		this->u2.setSize( this->mesh->getDofs() );
		this->p.setSize( this->mesh->getDofs() );
		this->e.setSize( this->mesh->getDofs() );
		}

		template< typename AdvectionScheme,
		@@ -226,6 +228,24 @@ const typename tnlNavierStokesSolver< AdvectionScheme, DiffusionScheme, Boundary
		return this->p;
		}

		template< typename AdvectionScheme,
		typename DiffusionScheme,
		typename BoundaryConditions >
		typename tnlNavierStokesSolver< AdvectionScheme, DiffusionScheme, BoundaryConditions >::VectorType&
		tnlNavierStokesSolver< AdvectionScheme, DiffusionScheme, BoundaryConditions >::getEnergy()
		{
		return this->e;
		}

		template< typename AdvectionScheme,
		typename DiffusionScheme,
		typename BoundaryConditions >
		const typename tnlNavierStokesSolver< AdvectionScheme, DiffusionScheme, BoundaryConditions >::VectorType&
		tnlNavierStokesSolver< AdvectionScheme, DiffusionScheme, BoundaryConditions >::getEnergy() const
		{
		return this->e;
		}

		template< typename AdvectionScheme,
		typename DiffusionScheme,
		typename BoundaryConditions >
		@@ -260,7 +280,8 @@ void tnlNavierStokesSolver< AdvectionScheme,
		DiffusionScheme,
		BoundaryConditions > :: updatePhysicalQuantities( const Vector& dofs_rho,
		const Vector& dofs_rho_u1,
		const Vector& dofs_rho_u2 )
		const Vector& dofs_rho_u2,
		const Vector& dofs_e )
		{
		if( DeviceType :: getDevice() == tnlHostDevice )
		{
		@@ -275,6 +296,7 @@ void tnlNavierStokesSolver< AdvectionScheme,
		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;
		this->e[ c ] = dofs_e[ c ];
		}
		}
		}
		@@ -295,25 +317,28 @@ void tnlNavierStokesSolver< AdvectionScheme,
		tnlAssert( this->u2Viscosity, );
		tnlAssert( this->boundaryConditions, );

		tnlSharedVector< RealType, DeviceType, IndexType > dofs_rho, dofs_rho_u1, dofs_rho_u2,
		rho_t, rho_u1_t, rho_u2_t;
		tnlSharedVector< RealType, DeviceType, IndexType > dofs_rho, dofs_rho_u1, dofs_rho_u2, dofs_e,
		rho_t, rho_u1_t, rho_u2_t, e_t;

		const IndexType& dofs = this->mesh->getDofs();
		dofs_rho. bind( & u. getData()[ 0 ], dofs );
		dofs_rho_u1. bind( & u. getData()[ dofs ], dofs );
		dofs_rho_u2. bind( & u. getData()[ 2 * dofs ], dofs );
		dofs_e. bind( & u. getData()[ 3 * dofs ], dofs );

		this->advection->setRho( dofs_rho );
		this->advection->setRhoU1( dofs_rho_u1 );
		this->advection->setRhoU2( dofs_rho_u2 );
		this->advection->setE( dofs_e );

		rho_t.bind( & fu. getData()[ 0 ], dofs );
		rho_u1_t.bind( & fu. getData()[ dofs ], dofs );
		rho_u2_t.bind( & fu. getData()[ 2 * dofs ], dofs );
		e_t.bind( & fu. getData()[ 3 * dofs ], dofs );

		updatePhysicalQuantities( dofs_rho, dofs_rho_u1, dofs_rho_u2 );
		updatePhysicalQuantities( dofs_rho, dofs_rho_u1, dofs_rho_u2, dofs_e );

		this->boundaryConditions->apply( time, rho, u1, u2 );
		this->boundaryConditions->apply( time, rho, u1, u2, e );

		const IndexType& xSize = this->mesh->getDimensions().x();
		const IndexType& ySize = this->mesh->getDimensions().y();
		@@ -330,9 +355,12 @@ void tnlNavierStokesSolver< AdvectionScheme,
		dofs_rho[ c1 ] = this->rho[ c1 ];
		dofs_rho_u1[ c1 ] = this->rho[ c1 ] * this->u1[ c1 ];
		dofs_rho_u2[ c1 ] = this->rho[ c1 ] * this->u2[ c1 ];
		dofs_e[ c1 ] = this->e[ c1 ];

		dofs_rho[ c3 ] = this->rho[ c3 ];
		dofs_rho_u1[ c3 ] = this->rho[ c3 ] * this->u1[ c3 ];
		dofs_rho_u2[ c3 ] = this->rho[ c3 ] * this->u2[ c3 ];
		dofs_e[ c3 ] = this->e[ c3 ];
		}
		for( IndexType j = 0; j < ySize; j ++ )
		{
		@@ -344,9 +372,12 @@ void tnlNavierStokesSolver< AdvectionScheme,
		dofs_rho[ c1 ] = this->rho[ c1 ];
		dofs_rho_u1[ c1 ] = this->rho[ c1 ] * this->u1[ c1 ];
		dofs_rho_u2[ c1 ] = this->rho[ c1 ] * this->u2[ c1 ];
		dofs_e[ c1 ] = this->e[ c1 ];

		dofs_rho[ c3 ] = this->rho[ c3 ];
		dofs_rho_u1[ c3 ] = this->rho[ c3 ] * this->u1[ c3 ];
		dofs_rho_u2[ c3 ] = this->rho[ c3 ] * this->u2[ c3 ];
		dofs_e[ c3 ] = this->e[ c3 ];
		}
		}

		@@ -369,6 +400,7 @@ void tnlNavierStokesSolver< AdvectionScheme,
		rho_t[ c ],
		rho_u1_t[ c ],
		rho_u2_t[ c ],
		e_t[ c ],
		tau );

		//rho_u1_t[ c ] += ;
		@@ -395,13 +427,14 @@ bool tnlNavierStokesSolver< AdvectionScheme,
		BoundaryConditions >::writePhysicalVariables( const RealType& t,
		const IndexType step )
		{
		tnlSharedVector< RealType, DeviceType, IndexType > dofs_rho, dofs_rho_u1, dofs_rho_u2;
		tnlSharedVector< RealType, DeviceType, IndexType > dofs_rho, dofs_rho_u1, dofs_rho_u2, dofs_e;
		const IndexType& dofs = mesh->getDofs();
		dofs_rho. bind( & dofVector.getData()[ 0 ], dofs );
		dofs_rho_u1. bind( & dofVector.getData()[ dofs ], dofs );
		dofs_rho_u2. bind( & dofVector.getData()[ 2 * dofs ], dofs );
		dofs_e. bind( & dofVector.getData()[ 3 * dofs ], dofs );

		this->updatePhysicalQuantities( dofs_rho, dofs_rho_u1, dofs_rho_u2 );
		this->updatePhysicalQuantities( dofs_rho, dofs_rho_u1, dofs_rho_u2, dofs_e );
		tnlVector< tnlTuple< 2, RealType >, DeviceType, IndexType > u;
		u. setLike( u1 );
		tnlString fileName;
		@@ -417,6 +450,9 @@ bool tnlNavierStokesSolver< AdvectionScheme,
		FileNameBaseNumberEnding( "p-", step, 5, ".tnl", fileName );
		if( ! this->p.save( fileName ) )
		return false;
		FileNameBaseNumberEnding( "e-", step, 5, ".tnl", fileName );
		if( ! this->e.save( fileName ) )
		return false;
		return true;
		}

src/solvers/cfd/navier-stokes/tnlNavierStokesSolver.h

+8 −2

Original line number	Diff line number	Diff line
		@@ -83,6 +83,11 @@ class tnlNavierStokesSolver

		const VectorType& getPressure() const;

		VectorType& getEnergy();

		const VectorType& getEnergy() const;


		IndexType getDofs() const;

		void bindDofVector( RealType* );
		@@ -92,7 +97,8 @@ class tnlNavierStokesSolver
		template< typename Vector >
		void updatePhysicalQuantities( const Vector& rho,
		const Vector& rho_u1,
		const Vector& rho_u2 );
		const Vector& rho_u2,
		const Vector& e );

		template< typename SolverVectorType >
		void getExplicitRhs( const RealType& time,
		@@ -117,7 +123,7 @@ class tnlNavierStokesSolver

		MeshType* mesh;

		VectorType rho, u1, u2, p;
		VectorType rho, u1, u2, p, e;

		RealType mu, gravity, R, T;