Implementing a non-linear diffusion to the Navier-Stokes solver. (413e639a) · Commits · TNL / tnl-dev

src/implementation/schemes/diffusion/tnlLinearDiffusion_impl.h

+80 −9

Original line number	Diff line number	Diff line
		@@ -267,7 +267,10 @@ void tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeomet
		}

		template< typename Real, typename Device, typename Index >
		Real tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry > > :: getDiffusion( const Index& c ) const
		Real tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry > > :: getDiffusion( const Index& c,
		const RealType t_11,
		const RealType t_22,
		const RealType t_12 ) const
		{
		tnlAssert( this -> mesh, cerr << "No mesh was set in tnlLinearDiffusion. Use the bindMesh method." );

		@@ -279,9 +282,77 @@ Real tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeomet
		const Index n = mesh->getElementNeighbour( c, 0, 1 );
		const Index s = mesh->getElementNeighbour( c, 0, -1 );

		return ( f[ e ] - 2.0 * f[ c ] + f[ w ] ) / ( hx * hx ) +
		( f[ n ] - 2.0 * f[ c ] + f[ s ] ) / ( hy * hy );
		RealType F_xx( 0.0 ), F_yy( 0.0 ), F_xy( 0.0 );
		if( t_11 != 0.0 )
		F_xx = ( f[ e ] - 2.0 * f[ c ] + f[ w ] ) / ( hx * hx );
		if( t_22 != 0.0 )
		F_yy = ( f[ n ] - 2.0 * f[ c ] + f[ s ] ) / ( hy * hy );
		if( t_12 != 0.0 )
		{
		const Index en = mesh->getElementNeighbour( c, 1, 1 );
		const Index wn = mesh->getElementNeighbour( c, -1, 1 );
		const Index es = mesh->getElementNeighbour( c, 1, -1 );
		const Index ws = mesh->getElementNeighbour( c, -1, -1 );

		const RealType F_en = 0.25 * ( f[ c ] + f[ e ] + f[ n ] + f[ en ] );
		const RealType F_wn = 0.25 * ( f[ c ] + f[ w ] + f[ n ] + f[ wn ] );
		const RealType F_es = 0.25 * ( f[ c ] + f[ e ] + f[ s ] + f[ es ] );
		const RealType F_ws = 0.25 * ( f[ c ] + f[ w ] + f[ s ] + f[ ws ] );

		const RealType F_y_e = ( F_en - F_es ) / hy;
		const RealType F_y_w = ( F_wn - F_ws ) / hy;
		F_xy = ( F_y_e - F_y_w ) / hx;
		}
		return t_11 * F_xx + t_22 * F_yy + t_12 * F_xy;
		}

		template< typename Real, typename Device, typename Index >
		template< typename Vector >
		Real tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeometry > >::getDiffusion( const IndexType& c,
		const Vector& v_11,
		const Vector& v_22,
		const Vector& v_12,
		RealType t_11,
		RealType t_22,
		RealType t_12 ) const
		{
		tnlAssert( this -> mesh, cerr << "No mesh was set in tnlLinearDiffusion. Use the bindMesh method." );

		const Real hx = mesh->getParametricStep().x();
		const Real hy = mesh->getParametricStep().y();

		const Index e = mesh->getElementNeighbour( c, 1, 0 );
		const Index w = mesh->getElementNeighbour( c, -1, 0 );
		const Index n = mesh->getElementNeighbour( c, 0, 1 );
		const Index s = mesh->getElementNeighbour( c, 0, -1 );

		const RealType F_x_e = ( f[ e ] - f[ c ] ) / hx;
		const RealType F_x_w = ( f[ c ] - f[ w ] ) / hx;
		const RealType F_y_n = ( f[ n ] - f[ c ] ) / hy;
		const RealType F_y_s = ( f[ c ] - f[ s ] ) / hy;

		RealType F_xx( 0.0 ), F_yy( 0.0 ), F_xy( 0.0 );
		if( t_11 != 0.0 )
		F_xx = ( v_11[ e ] * f[ e ] - v_11[ w ] * f[ w ] ) / hx;
		if( t_22 != 0.0 )
		F_yy = ( v_22[ n ] * f[ n ] - v_22[ s ] * f[ s ] ) / hy;
		if( t_12 != 0.0 )
		{
		const Index en = mesh->getElementNeighbour( c, 1, 1 );
		const Index wn = mesh->getElementNeighbour( c, -1, 1 );
		const Index es = mesh->getElementNeighbour( c, 1, -1 );
		const Index ws = mesh->getElementNeighbour( c, -1, -1 );

		const RealType F_en = 0.25 * ( f[ c ] + f[ e ] + f[ n ] + f[ en ] );
		const RealType F_wn = 0.25 * ( f[ c ] + f[ w ] + f[ n ] + f[ wn ] );
		const RealType F_es = 0.25 * ( f[ c ] + f[ e ] + f[ s ] + f[ es ] );
		const RealType F_ws = 0.25 * ( f[ c ] + f[ w ] + f[ s ] + f[ ws ] );

		const RealType F_y_e = ( F_en - F_es ) / hy;
		const RealType F_y_w = ( F_wn - F_ws ) / hy;
		const RealType F_xy = ( v_12[ e ] * F_y_e - v_12[ w ] * F_y_w ) / hx;
		}
		return t_11 * F_xx + t_12 * F_xy + t_22 * F_yy;
		}

		#endif

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

+26 −27

Original line number	Diff line number	Diff line
		@@ -316,9 +316,9 @@ void tnlNavierStokesSolver< AdvectionScheme,
		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;
		//this->p[ c ] = ( this->gamma - 1.0 ) *
		// ( dofs_e[ c ] - 0.5 * this->rho[ c ] * ( this->u1[ c ] * this->u1[ c ] + this->u2[ c ] * this->u2[ c ] ) );
		//this->p[ c ] = dofs_rho[ c ] * this -> R * this -> T;
		this->p[ c ] = ( this->gamma - 1.0 ) *
		( dofs_e[ c ] - 0.5 * this->rho[ c ] * ( this->u1[ c ] * this->u1[ c ] + this->u2[ c ] * this->u2[ c ] ) );
		this->temperature[ c ] = this->p[ c ] / ( this->rho[ c ] * this->R );
		}
		}
		@@ -452,32 +452,31 @@ void tnlNavierStokesSolver< AdvectionScheme,
		/***
		* Add the viscosity term
		*/
		rho_u1_t[ c ] += this->mu*u1Viscosity->getDiffusion( c );;
		rho_u2_t[ c ] += this->mu*u2Viscosity->getDiffusion( c );

		IndexType e = this->mesh->getElementIndex( i+1, j );
		IndexType w = this->mesh->getElementIndex( i-1, j );
		IndexType n = this->mesh->getElementIndex( i, j+1 );
		IndexType s = this->mesh->getElementIndex( i, j-1 );
		RealType u1_e = 0.5*( this->u1[ c ] + this->u1[ e ] );
		RealType u1_w = 0.5*( this->u1[ c ] + this->u1[ w ] );
		RealType u2_n = 0.5*( this->u2[ c ] + this->u2[ n ] );
		RealType u2_s = 0.5*( this->u2[ c ] + this->u2[ s ] );
		RealType hx = mesh->getParametricStep().x();
		RealType hy = mesh->getParametricStep().y();
		RealType u1_x_e = ( this->u1[ e ] - this->u1[ c ] ) / hx;
		RealType u1_x_w = ( this->u1[ c ] - this->u1[ w ] ) / hx;
		RealType u2_y_n = ( this->u2[ n ] - this->u2[ c ] ) / hy;
		RealType u2_y_s = ( this->u2[ c ] - this->u2[ s ] ) / hy;


		e_t[ c ] += this->mu( ( u1_e u1_x_e - u1_w * u1_x_w ) / hx +
		( u2_n * u2_y_n - u2_s * u2_y_s ) / hy );
		e_t[ c ] = 0.0;
		}

		writeExplicitRhs( time, -1 );
		getchar();
		rho_u1_t[ c ] += this->mu*( u1Viscosity->getDiffusion( c, 4.0/3.0, 1.0, 0.0 ) +
		u2Viscosity->getDiffusion( c, 0.0, 0.0, 1.0/3.0 ) );
		rho_u2_t[ c ] += this->mu*( u2Viscosity->getDiffusion( c, 1.0, 4.0/3.0, 0.0 ) +
		u1Viscosity->getDiffusion( c, 0.0, 0.0, 1.0/3.0 ) );


		RealType k = 2.495*pow( 400.0, 1.5 ) / ( 400.0 + 194.0 );
		e_t[ c ] += this->mu*( u1Viscosity->getDiffusion( c,
		this->u1, this->u1, this->u2,
		4.0/3.0, 1.0, -2.0/3.0 ) +
		u1Viscosity->getDiffusion( c,
		this->u1, this->u1, this->u2,
		0.0, 0.0, 1.0 ) +
		u2Viscosity->getDiffusion( c,
		this->u2, this->u2, this->u1,
		1.0, 4.0/3.0, -2.0/3.0 ) +
		u2Viscosity->getDiffusion( c,
		this->u2, this->u2, this->u1,
		0.0, 0.0, 1.0 ) +
		k * temperatureViscosity->getDiffusion( c, 1.0, 1.0, 0.0 ) );
		//e_t[ c ] = 0.0;
		}

		//writeExplicitRhs( time, -1 );
		//getchar();
		}

		template< typename AdvectionScheme,

src/schemes/diffusion/tnlLinearDiffusion.h

+14 −1

Original line number	Diff line number	Diff line
		@@ -77,7 +77,20 @@ class tnlLinearDiffusion< tnlGrid< 2, Real, Device, Index, tnlIdenticalGridGeome
		template< typename Vector >
		void setFunction( Vector& f ); // TODO: add const

		RealType getDiffusion( const Index& i ) const;
		RealType getDiffusion( const IndexType& i,
		const RealType t_11 = 1.0,
		const RealType t_22 = 1.0,
		const RealType t_12 = 0.0 ) const;

		template< typename Vector >
		RealType getDiffusion( const IndexType& i,
		const Vector& v_11,
		const Vector& v_22,
		const Vector& v_12,
		RealType t_11,
		RealType t_22,
		RealType t_12 ) const;

		protected:

		// TODO: change to ConstSharedVector