Fixed FSM interface initiation in CUDA. (6d1a8d98) · Commits · TNL / tnl-dev

src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h

+18 −12

Original line number	Diff line number	Diff line
		@@ -31,10 +31,12 @@ class tnlDirectEikonalMethodsBase< Meshes::Grid< 1, Real, Device, Index > >
		typedef Index IndexType;
		typedef Functions::MeshFunction< MeshType > MeshFunctionType;
		typedef Functions::MeshFunction< MeshType, 1, bool > InterfaceMapType;
		using MeshFunctionPointer = SharedPointer< MeshFunctionType >;
		using InterfaceMapPointer = SharedPointer< InterfaceMapType >;

		void initInterface( const MeshFunctionType& input,
		MeshFunctionType& output,
		InterfaceMapType& interfaceMap );
		void initInterface( const MeshFunctionPointer& input,
		MeshFunctionPointer& output,
		InterfaceMapPointer& interfaceMap );

		template< typename MeshEntity >
		__cuda_callable__ void updateCell( MeshFunctionType& u,
		@@ -55,10 +57,12 @@ class tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > >
		typedef Index IndexType;
		typedef Functions::MeshFunction< MeshType > MeshFunctionType;
		typedef Functions::MeshFunction< MeshType, 2, bool > InterfaceMapType;
		using MeshFunctionPointer = SharedPointer< MeshFunctionType >;
		using InterfaceMapPointer = SharedPointer< InterfaceMapType >;

		void initInterface( const MeshFunctionType& input,
		MeshFunctionType& output,
		InterfaceMapType& interfaceMap );
		void initInterface( const MeshFunctionPointer& input,
		MeshFunctionPointer& output,
		InterfaceMapPointer& interfaceMap );

		template< typename MeshEntity >
		__cuda_callable__ void updateCell( MeshFunctionType& u,
		@@ -81,10 +85,12 @@ class tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > >
		typedef Index IndexType;
		typedef Functions::MeshFunction< MeshType > MeshFunctionType;
		typedef Functions::MeshFunction< MeshType, 3, bool > InterfaceMapType;
		using MeshFunctionPointer = SharedPointer< MeshFunctionType >;
		using InterfaceMapPointer = SharedPointer< InterfaceMapType >;

		void initInterface( const MeshFunctionType& input,
		MeshFunctionType& output,
		InterfaceMapType& interfaceMap );
		void initInterface( const MeshFunctionPointer& input,
		MeshFunctionPointer& output,
		InterfaceMapPointer& interfaceMap );

		template< typename MeshEntity >
		__cuda_callable__ void updateCell( MeshFunctionType& u,
		@@ -109,12 +115,12 @@ template < typename Real, typename Device, typename Index >
		__global__ void CudaUpdateCellCaller( Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index >, 2, bool >& interfaceMap,
		Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index > >& aux );

		/*template < typename Real, typename Device, typename Index >
		template < typename Real, typename Device, typename Index >
		__global__ void CudaInitCaller( const Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index > >& input,
		Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index > >& output,
		Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index >, 2, bool >& interfaceMap );*/
		Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index >, 2, bool >& interfaceMap );

		__global__ void CudaInitCaller( const Functions::MeshFunction< Meshes::Grid< 2, double, TNL::Devices::Cuda, int > >& input );
		//__global__ void CudaInitCaller( const Functions::MeshFunction< Meshes::Grid< 2, double, TNL::Devices::Cuda, int > >& input );
		#endif

		#include "tnlDirectEikonalMethodsBase_impl.h"

src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h

+30 −19

Original line number	Diff line number	Diff line
		@@ -16,12 +16,15 @@ template< typename Real,
		typename Index >
		void
		tnlDirectEikonalMethodsBase< Meshes::Grid< 1, Real, Device, Index > >::
		initInterface( const MeshFunctionType& input,
		MeshFunctionType& output,
		InterfaceMapType& interfaceMap )
		initInterface( const MeshFunctionPointer& _input,
		MeshFunctionPointer& _output,
		InterfaceMapPointer& _interfaceMap )
		{
		const MeshType& mesh = input.getMesh();
		const MeshType& mesh = _input->getMesh();
		typedef typename MeshType::Cell Cell;
		const MeshFunctionType& input = _input.getData();
		MeshFunctionType& output = _output.modifyData();
		InterfaceMapType& interfaceMap = _interfaceMap.modifyData();
		Cell cell( mesh );
		for( cell.getCoordinates().x() = 1;
		cell.getCoordinates().x() < mesh.getDimensions().x() - 1;
		@@ -74,9 +77,9 @@ template< typename Real,
		typename Index >
		void
		tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > >::
		initInterface( const MeshFunctionType& input,
		MeshFunctionType& output,
		InterfaceMapType& interfaceMap )
		initInterface( const MeshFunctionPointer& _input,
		MeshFunctionPointer& _output,
		InterfaceMapPointer& _interfaceMap )
		{
		/*
		doplnit přepočty pro cudu:
		@@ -88,7 +91,7 @@ initInterface( const MeshFunctionType& input,
		if( std::is_same< Device, Devices::Cuda >::value )
		{
		#ifdef HAVE_CUDA
		const MeshType& mesh = input.getMesh();
		const MeshType& mesh = _input->getMesh();

		const int cudaBlockSize( 16 );
		int numBlocksX = Devices::Cuda::getNumberOfBlocks( mesh.getDimensions().x(), cudaBlockSize );
		@@ -96,14 +99,19 @@ initInterface( const MeshFunctionType& input,
		dim3 blockSize( cudaBlockSize, cudaBlockSize );
		dim3 gridSize( numBlocksX, numBlocksY );
		Devices::Cuda::synchronizeDevice();
		//CudaInitCaller< Real, Device, Index ><<< gridSize, blockSize >>>( input, output, interfaceMap );
		CudaInitCaller<<< gridSize, blockSize >>>( input );
		CudaInitCaller<<< gridSize, blockSize >>>( _input.template getData< Device >(),
		_output.template modifyData< Device >(),
		_interfaceMap.template modifyData< Device >() );
		//CudaInitCaller<<< gridSize, blockSize >>>( input );
		cudaDeviceSynchronize();
		TNL_CHECK_CUDA_DEVICE;
		#endif
		}
		if( std::is_same< Device, Devices::Host >::value )
		{
		const MeshFunctionType& input = _input.getData();
		MeshFunctionType& output = _output.modifyData();
		InterfaceMapType& interfaceMap = _interfaceMap.modifyData();
		const MeshType& mesh = input.getMesh();
		typedef typename MeshType::Cell Cell;
		Cell cell( mesh );
		@@ -289,10 +297,13 @@ template< typename Real,
		typename Index >
		void
		tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > >::
		initInterface( const MeshFunctionType& input,
		MeshFunctionType& output,
		InterfaceMapType& interfaceMap )
		initInterface( const MeshFunctionPointer& _input,
		MeshFunctionPointer& _output,
		InterfaceMapPointer& _interfaceMap )
		{
		const MeshFunctionType& input = _input.getData();
		MeshFunctionType& output = _output.modifyData();
		InterfaceMapType& interfaceMap = _interfaceMap.modifyData();
		const MeshType& mesh = input.getMesh();
		typedef typename MeshType::Cell Cell;
		Cell cell( mesh );
		@@ -602,16 +613,16 @@ __cuda_callable__ void sortMinims( T1 pom[])
		}
		}

		/*template < typename Real, typename Device, typename Index >
		template < typename Real, typename Device, typename Index >
		__global__ void CudaInitCaller( const Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index > >& input,
		Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index > >& output,
		Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index >, 2, bool >& interfaceMap )
		{
		int i = threadIdx.x + blockDim.x*blockIdx.x;
		int j = blockDim.y*blockIdx.y + threadIdx.y;
		const Meshes::Grid< 2, Real, Device, Index >& mesh = input.getMesh();
		const Meshes::Grid< 2, Real, Device, Index >& mesh = input.template getMesh< Devices::Cuda >();

		//if( i < mesh.getDimensions().x() && j < mesh.getDimensions().y() )
		if( i < mesh.getDimensions().x() && j < mesh.getDimensions().y() )
		{
		typedef typename Meshes::Grid< 2, Real, Device, Index >::Cell Cell;
		Cell cell( mesh );
		@@ -671,13 +682,13 @@ __global__ void CudaInitCaller( const Functions::MeshFunction< Meshes::Grid< 2,
		}
		}
		}
		}*/
		}


		__global__ void CudaInitCaller( const Functions::MeshFunction< Meshes::Grid< 2, double, TNL::Devices::Cuda, int > >& input )
		/*__global__ void CudaInitCaller( const Functions::MeshFunction< Meshes::Grid< 2, double, TNL::Devices::Cuda, int > >& input )
		{
		int i = threadIdx.x + blockDim.x*blockIdx.x;
		int j = blockDim.y*blockIdx.y + threadIdx.y;
		//const Meshes::Grid< 2, double, TNL::Devices::Cuda, int >& mesh = input.getMesh();

		}
		}*/

src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem.h

+5 −3

Original line number	Diff line number	Diff line
		@@ -36,6 +36,8 @@ class tnlDirectEikonalProblem
		typedef Functions::MeshFunction< Mesh > MeshFunctionType;
		typedef Problems::PDEProblem< Mesh, RealType, DeviceType, IndexType > BaseType;
		using AnisotropyType = Anisotropy;
		using AnisotropyPointer = SharedPointer< AnisotropyType, DeviceType >;
		using MeshFunctionPointer = SharedPointer< MeshFunctionType >;

		using typename BaseType::MeshType;
		using typename BaseType::DofVectorType;
		@@ -76,11 +78,11 @@ class tnlDirectEikonalProblem

		protected:

		MeshFunctionType u;
		MeshFunctionPointer u;

		MeshFunctionType initialData;
		MeshFunctionPointer initialData;

		AnisotropyType anisotropy;
		AnisotropyPointer anisotropy;

		};

src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem_impl.h

+3 −3

Original line number	Diff line number	Diff line
		@@ -95,7 +95,7 @@ tnlDirectEikonalProblem< Mesh, Anisotropy, Real, Index >::
		bindDofs( const MeshPointer& mesh,
		const DofVectorPointer& dofs )
		{
		this->u.bind( mesh, dofs );
		this->u->bind( mesh, dofs );
		}

		template< typename Mesh,
		@@ -110,8 +110,8 @@ setInitialCondition( const Config::ParameterContainer& parameters,
		MeshDependentDataPointer& meshdependentData )
		{
		String inputFile = parameters.getParameter< String >( "input-file" );
		this->initialData.setMesh( mesh );
		if( !this->initialData.boundLoad( inputFile ) )
		this->initialData->setMesh( mesh );
		if( !this->initialData->boundLoad( inputFile ) )
		return false;
		return true;
		}

src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod.h

+19 −11

Original line number	Diff line number	Diff line
		@@ -34,14 +34,18 @@ class FastSweepingMethod< Meshes::Grid< 1, Real, Device, Index >, Anisotropy >

		typedef Meshes::Grid< 1, Real, Device, Index > MeshType;
		typedef Real RealType;
		typedef TNL::Devices::Host DeviceType;
		typedef Device DeviceType;
		typedef Index IndexType;
		typedef Anisotropy AnisotropyType;
		typedef tnlDirectEikonalMethodsBase< Meshes::Grid< 1, Real, Device, Index > > BaseType;
		using MeshPointer = SharedPointer< MeshType >;
		using AnisotropyPointer = SharedPointer< AnisotropyType, DeviceType >;

		using typename BaseType::InterfaceMapType;
		using typename BaseType::MeshFunctionType;
		using typename BaseType::InterfaceMapPointer;
		using typename BaseType::MeshFunctionPointer;


		FastSweepingMethod();

		@@ -50,8 +54,8 @@ class FastSweepingMethod< Meshes::Grid< 1, Real, Device, Index >, Anisotropy >
		void setMaxIterations( const IndexType& maxIterations );

		void solve( const MeshPointer& mesh,
		const AnisotropyType& anisotropy,
		MeshFunctionType& u );
		const AnisotropyPointer& anisotropy,
		MeshFunctionPointer& u );


		protected:
		@@ -72,14 +76,17 @@ class FastSweepingMethod< Meshes::Grid< 2, Real, Device, Index >, Anisotropy >

		typedef Meshes::Grid< 2, Real, Device, Index > MeshType;
		typedef Real RealType;
		typedef TNL::Devices::Host DeviceType;
		typedef Device DeviceType;
		typedef Index IndexType;
		typedef Anisotropy AnisotropyType;
		typedef tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > > BaseType;
		using MeshPointer = SharedPointer< MeshType >;
		using AnisotropyPointer = SharedPointer< AnisotropyType, DeviceType >;

		using typename BaseType::InterfaceMapType;
		using typename BaseType::MeshFunctionType;
		using typename BaseType::InterfaceMapPointer;
		using typename BaseType::MeshFunctionPointer;

		FastSweepingMethod();

		@@ -88,8 +95,8 @@ class FastSweepingMethod< Meshes::Grid< 2, Real, Device, Index >, Anisotropy >
		void setMaxIterations( const IndexType& maxIterations );

		void solve( const MeshPointer& mesh,
		const AnisotropyType& anisotropy,
		MeshFunctionType& u );
		const AnisotropyPointer& anisotropy,
		MeshFunctionPointer& u );


		protected:
		@@ -110,14 +117,17 @@ class FastSweepingMethod< Meshes::Grid< 3, Real, Device, Index >, Anisotropy >

		typedef Meshes::Grid< 3, Real, Device, Index > MeshType;
		typedef Real RealType;
		typedef TNL::Devices::Host DeviceType;
		typedef Device DeviceType;
		typedef Index IndexType;
		typedef Anisotropy AnisotropyType;
		typedef tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > > BaseType;
		using MeshPointer = SharedPointer< MeshType >;
		using AnisotropyPointer = SharedPointer< AnisotropyType, DeviceType >;

		using typename BaseType::InterfaceMapType;
		using typename BaseType::MeshFunctionType;
		using typename BaseType::InterfaceMapPointer;
		using typename BaseType::MeshFunctionPointer;

		FastSweepingMethod();

		@@ -126,8 +136,8 @@ class FastSweepingMethod< Meshes::Grid< 3, Real, Device, Index >, Anisotropy >
		void setMaxIterations( const IndexType& maxIterations );

		void solve( const MeshPointer& mesh,
		const AnisotropyType& anisotropy,
		MeshFunctionType& u );
		const AnisotropyPointer& anisotropy,
		MeshFunctionPointer& u );


		protected:
		@@ -135,8 +145,6 @@ class FastSweepingMethod< Meshes::Grid< 3, Real, Device, Index >, Anisotropy >
		const IndexType maxIterations;
		};



		#include "tnlFastSweepingMethod1D_impl.h"
		#include "tnlFastSweepingMethod2D_impl.h"
		#include "tnlFastSweepingMethod3D_impl.h"