Working in the Navier-Stokes solver. (1ea595f6) · Commits · TNL / tnl-dev

examples/navier-stokes/Makefile

+4 −4

Original line number	Diff line number	Diff line
		@@ -9,14 +9,14 @@ 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) -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) -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
		#LD_FLAGS = -L$(TNL_INSTALL_DIR) -ltnl-dbg-0.1 -lgomp
		#LD_FLAGS = -L$(TNL_INSTALL_DIR) -ltnl-0.1 -lgomp
		LD_FLAGS = -L$(TNL_INSTALL_DIR) -ltnl-dbg-0.1 -lgomp

		SOURCES = main.cpp
		HEADERS = navierStokesSetter.h \

examples/navier-stokes/make-png-vectors-from-gnuplot

0 → 100644

+75 −0

Original line number	Diff line number	Diff line
		#!/bin/bash

		# $1 - base file name
		# $2 - xrange
		# $3 - yrange
		# $4 - cbrange
		# $5 - vector field skipping along x
		# $6 - vector field skipping along y

		function processFile()
		{
		file=${1}

		gnuplotcommand="
		set terminal png giant size 1280,1280 crop;
		set output '`basename $file ".gplt"`.png';
		set pm3d map;
		set palette defined(0.0 0.5 1.0 0, 0.02 \"light-goldenrod\", 0.04 \"yellow\", 0.08 \"red\", 0.4 \"light-blue\", 1.0 \"blue\");
		unset key;
		set size ratio -1;
		set pointsize 0.4;"

		if ! test x$2 = x;
		then
		gnuplotcommand="${gnuplotcommand} set xrange [0:$2];"
		fi
		if ! test x$3 = x;
		then
		gnuplotcommand="${gnuplotcommand} set yrange [0:$3];"
		fi
		if ! test x$4 = x;
		then
		gnuplotcommand="${gnuplotcommand} set cbrange [0:$4];"
		fi

		gnuplotcommand="${gnuplotcommand} splot '$file' using 1:2:(sqrt(\$32 + \$42)) w pm3d title '${1}';"
		echo ${gnuplotcommand} \| gnuplot


		gnuplotcommand="
		set terminal png giant size 1280,1280 crop;
		set output 'vec-`basename $file ".gplt"`.png';
		set palette defined(0.0 0.5 1.0 0, 0.02 \"light-goldenrod\", 0.04 \"yellow\", 0.08 \"red\", 0.4 \"light-blue\", 1.0 \"blue\");
		unset key;
		set size ratio -1;
		set pointsize 0.4;"

		if ! test x$2 = x;
		then
		gnuplotcommand="${gnuplotcommand} set xrange [0:$2];"
		fi
		if ! test x$3 = x;
		then
		gnuplotcommand="${gnuplotcommand} set yrange [0:$3];"
		fi
		if ! test x$4 = x;
		then
		gnuplotcommand="${gnuplotcommand} set cbrange [0:$4];"
		fi

		gnuplotcommand="${gnuplotcommand} plot '$file' every ${5}:${6} using 1:2:3:4:(sqrt($32+$42)) with vectors linecolor palette z title '${1}';"
		echo ${gnuplotcommand} \| gnuplot
		}

		for file in ${1}*.gplt
		do
		png_file="`basename $file ".gplt"`.png"
		if test -e ${png_file};
		then
		echo -ne "Skipping: ${png_file} \r"
		else
		echo -ne "Creating: ${png_file} \r"
		processFile ${file} ${2} ${3} ${4} ${5} ${6}
		fi
		done

examples/navier-stokes/navierStokesSolver.h

+5 −5

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

		#include "navierStokesSolverMonitor.h"
		#include "navierStokesBoundaryConditions.h"
		@@ -102,11 +102,11 @@ class navierStokesSolver

		DofVectorType dofVector, rhsDofVector;

		RealType mu, R, T, p_0, gravity;
		RealType p_0, gravity;

		EulerScheme eulerScheme;

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

examples/navier-stokes/navierStokesSolver_impl.h

+15 −183

Original line number	Diff line number	Diff line
		@@ -47,9 +47,7 @@ __device__ void computeVelocityFieldCuda( const Index size,

		template< typename Mesh, typename EulerScheme >
		navierStokesSolver< Mesh, EulerScheme > :: navierStokesSolver()
		: R( 0.0 ),
		T( 0.0 ),
		p_0( 0.0 )
		: p_0( 0.0 )
		{

		this -> mesh. setName( "navier-stokes-mesh" );
		@@ -153,6 +151,7 @@ bool navierStokesSolver< Mesh, EulerScheme > :: init( const tnlParameterContaine
		/****
		* Set-up numerical scheme
		*/
		navierStokesScheme.setMesh( this->mesh );
		pressureGradient.setFunction( navierStokesScheme.getPressure() );
		pressureGradient.bindMesh( this -> mesh );
		this->eulerScheme. bindMesh( this -> mesh );
		@@ -162,7 +161,11 @@ bool navierStokesSolver< Mesh, EulerScheme > :: init( const tnlParameterContaine
		this->u2Viscosity. bindMesh( this -> mesh );
		this->u2Viscosity. setFunction( this -> navierStokesScheme.getU2() );
		navierStokesScheme.setAdvectionScheme( this->eulerScheme );
		navierStokesScheme.setMesh( this->mesh );
		navierStokesScheme.setDiffusionScheme( this->u1Viscosity,
		this->u2Viscosity );
		navierStokesScheme.setBoundaryConditions( this->boundaryConditions );
		navierStokesScheme.bindDofVector( dofVector.getData() );

		//navierStokesScheme.setDifusionScheme( this)
		return true;
		}
		@@ -176,7 +179,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 -> R * this -> T ) );
		dofs_rho. setValue( p_0 / ( this->navierStokesScheme.getR() *
		this->navierStokesScheme.getT() ) );
		rho_u1. setValue( 0.0 );
		rho_u2. setValue( 0.0 );

		@@ -192,7 +196,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 -> R * this -> T ) );
		dofs_rho. setElement( c, p_0 / ( this->navierStokesScheme.getR() *
		this->navierStokesScheme.getT() ) );
		rho_u1. setElement( c, 0.0 );
		rho_u2. setElement( c, 0.0 );

		@@ -211,193 +216,20 @@ bool navierStokesSolver< Mesh, EulerScheme > :: makeSnapshot( const RealType& t,
		const IndexType step )
		{
		cout << endl << "Writing output at time " << t << " step " << step << "." << endl;
		tnlSharedVector< RealType, DeviceType, IndexType > rho, rho_u1, rho_u2,
		rho_t, rho_u1_t, rho_u2_t;
		const IndexType& dofs = mesh. getDofs();
		rho. bind( & dofVector.getData()[ 0 ], dofs );
		rho_u1. bind( & dofVector.getData()[ dofs ], dofs );
		rho_u2. bind( & dofVector.getData()[ 2 * dofs ], dofs );
		rho_t. bind( & this->rhsDofVector.getData()[ 0 ], dofs );
		rho_u1_t. bind( & this->rhsDofVector.getData()[ dofs ], dofs );
		rho_u2_t. bind( & this->rhsDofVector.getData()[ 2 * dofs ], dofs );


		navierStokesScheme.updatePhysicalQuantities( rho, rho_u1, rho_u2 );
		tnlVector< RealType, DeviceType, IndexType > u;
		u. setLike( navierStokesScheme.getU1() );
		for( IndexType i = 0; i < navierStokesScheme.getU1().getSize(); 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 ) )
		return false;
		FileNameBaseNumberEnding( "u-2-", step, 5, ".tnl", fileName );
		if( ! this -> u2. save( fileName ) )
		return false;*/
		/*FileNameBaseNumberEnding( "p-", step, 5, ".tnl", fileName );
		if( ! p. save( fileName ) )
		return false;*/
		FileNameBaseNumberEnding( "u-", step, 5, ".tnl", fileName );
		if( ! u. save( fileName ) )
		return false;
		/*FileNameBaseNumberEnding( "rho-t-", step, 5, ".tnl", fileName );
		if( ! rho_t. save( fileName ) )
		return false;
		FileNameBaseNumberEnding( "rho-u1-t-", step, 5, ".tnl", fileName );
		if( ! rho_u1_t. save( fileName ) )
		if( !navierStokesScheme.writePhysicalVariables( t, step ) )
		return false;

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

		FileNameBaseNumberEnding( "rho-", step, 5, ".tnl", fileName );
		if( ! navierStokesScheme.getRho(). save( fileName ) )
		return false;
		FileNameBaseNumberEnding( "rho-u1-", step, 5, ".tnl", fileName );
		if( ! rho_u1. save( fileName ) )
		return false;

		FileNameBaseNumberEnding( "rho-u2-", step, 5, ".tnl", fileName );
		if( ! rho_u2. save( fileName ) )
		if( !navierStokesScheme.writeConservativeVariables( t, step ) )
		return false;
		return true;
		}

		/*template< typename Mesh, typename EulerScheme >
		template< typename Vector >
		void navierStokesSolver< Mesh, EulerScheme > :: updatePhysicalQuantities( const Vector& dofs_rho,
		const Vector& rho_u1,
		const Vector& rho_u2 )
		{
		if( DeviceType :: getDevice() == tnlHostDevice )
		{
		const IndexType& xSize = mesh. getDimensions(). x();
		const IndexType& ySize = mesh. getDimensions(). y();

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

		template< typename Mesh, typename EulerScheme >
		void navierStokesSolver< Mesh, EulerScheme > :: GetExplicitRHS( const RealType& time,
		const RealType& tau,
		DofVectorType& u,
		DofVectorType& fu )
		{
		tnlSharedVector< RealType, DeviceType, IndexType > dofs_rho, rho_u1, rho_u2,
		rho_t, rho_u1_t, rho_u2_t;

		const IndexType& dofs = mesh. getDofs();
		dofs_rho. bind( & u. getData()[ 0 ], dofs );
		rho_u1. bind( & u. getData()[ dofs ], dofs );
		rho_u2. bind( & u. getData()[ 2 * dofs ], dofs );

		eulerScheme. setRho( dofs_rho );
		eulerScheme. setRhoU1( rho_u1 );
		eulerScheme. setRhoU2( rho_u2 );

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

		navierStokesScheme.updatePhysicalQuantities( dofs_rho, rho_u1, rho_u2 );

		/****
		* Compute RHS
		*/
		const IndexType& xSize = mesh. getDimensions(). x();
		const IndexType& ySize = mesh. getDimensions(). y();
		/*const RealType hx = mesh. getParametricStep(). x();
		const RealType hy = mesh. getParametricStep(). y();*/

		if( DeviceType :: getDevice() == tnlHostDevice )
		{
		this->boundaryConditions.apply( time,
		this->navierStokesScheme.getRho(),
		this->navierStokesScheme.getU1(),
		this->navierStokesScheme.getU2() );
		for( IndexType i = 0; i < xSize; 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 );

		dofs_rho[ c1 ] = this->navierStokesScheme.getRho()[ c1 ];
		rho_u1[ c1 ] = this->navierStokesScheme.getRho()[ c1 ] * this->navierStokesScheme.getU1()[ c1 ];
		rho_u2[ c1 ] = this->navierStokesScheme.getRho()[ c1 ] * this->navierStokesScheme.getU2()[ c1 ];
		dofs_rho[ c3 ] = this->navierStokesScheme.getRho()[ c3 ];
		rho_u1[ c3 ] = this->navierStokesScheme.getRho()[ c3 ] * this->navierStokesScheme.getU1()[ c3 ];
		rho_u2[ c3 ] = this->navierStokesScheme.getRho()[ c3 ] * this->navierStokesScheme.getU2()[ c3 ];
		}
		for( IndexType j = 0; j < ySize; j ++ )
		{
		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 );

		dofs_rho[ c1 ] = this->navierStokesScheme.getRho()[ c1 ];
		rho_u1[ c1 ] = this->navierStokesScheme.getRho()[ c1 ] * this->navierStokesScheme.getU1()[ c1 ];
		rho_u2[ c1 ] = this->navierStokesScheme.getRho()[ c1 ] * this->navierStokesScheme.getU2()[ c1 ];
		dofs_rho[ c3 ] = this->navierStokesScheme.getRho()[ c3 ];
		rho_u1[ c3 ] = this->navierStokesScheme.getRho()[ c3 ] * this->navierStokesScheme.getU1()[ c3 ];
		rho_u2[ c3 ] = this->navierStokesScheme.getRho()[ c3 ] * this->navierStokesScheme.getU2()[ c3 ];

		}

		#ifdef HAVE_OPENMP
		#pragma omp parallel for
		#endif
		for( IndexType j = 0; j < ySize; j ++ )
		for( IndexType i = 0; i < xSize; i ++ )
		{
		IndexType c = mesh. getElementIndex( i, j );
		if( i == 0 \|\| j == 0 \|\|
		i == xSize - 1 \|\| j == ySize - 1 )
		{
		rho_t[ c ] = rho_u1_t[ c ] = rho_u2_t[ c ] = 0.0;
		continue;
		}

		eulerScheme. getExplicitRhs( c,
		rho_t[ c ],
		rho_u1_t[ c ],
		rho_u2_t[ c ],
		tau );

		//rho_u1_t[ c ] += ;
		//rho_u2_t[ c ] -= startUpCoefficient * this -> gravity * this -> rho[ c ];

		/***
		* Add the viscosity term
		*/
		rho_u1_t[ c ] += this->navierStokesScheme.getMu() * u1Viscosity. getDiffusion( c );
		rho_u2_t[ c ] += this->navierStokesScheme.getMu() * u2Viscosity. getDiffusion( c );

		}
		}

		/*rhsDofVector = fu;
		makeSnapshot( 0.0, 1 );
		getchar();*/

		navierStokesScheme.getExplicitRhs( time, tau, u, fu );
		}

		template< typename Mesh, typename EulerScheme >

examples/navier-stokes/share/examples/cavity

+3 −2

Original line number	Diff line number	Diff line
		@@ -39,13 +39,14 @@ do
		--verbose 2 \
		--device host

		tnl-view --mesh mesh.tnl --input-files u*tnl
		tnl-view --mesh mesh.tnl --check-output-file yes --input-files u*tnl

		make-png-from-gnuplot u 1 1 1.5
		make-png-vectors-from-gnuplot u 1 1 1.5 1 1

		rm *.gplt

		mencoder "mf://u-*.png" -mf fps=25 -o u.avi -ovc lavc -lavcopts vcodec=mpeg4
		mencoder "mf://vec-u-*.png" -mf fps=25 -o vec-u.avi -ovc lavc -lavcopts vcodec=mpeg4

		rm *.png