From 6d1a8d988ab4960715d100af5c6023f3b719207f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= <oberhuber.tomas@gmail.com> Date: Wed, 18 Apr 2018 21:19:42 +0200 Subject: [PATCH] Fixed FSM interface initiation in CUDA. --- .../tnlDirectEikonalMethodsBase.h | 30 +++++++----- .../tnlDirectEikonalMethodsBase_impl.h | 49 ++++++++++++------- .../hamilton-jacobi/tnlDirectEikonalProblem.h | 8 +-- .../tnlDirectEikonalProblem_impl.h | 6 +-- .../hamilton-jacobi/tnlFastSweepingMethod.h | 30 +++++++----- .../tnlFastSweepingMethod1D_impl.h | 14 +++--- .../tnlFastSweepingMethod2D_impl.h | 18 ++++--- .../tnlFastSweepingMethod3D_impl.h | 18 ++++--- 8 files changed, 102 insertions(+), 71 deletions(-) diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h index ed176f08a9..d8770b0e73 100644 --- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h +++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h @@ -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" diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h index 8444b6c378..37fd245128 100644 --- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h +++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h @@ -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(); -} +}*/ diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem.h index ca3b8368cc..66513d76f9 100644 --- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem.h +++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem.h @@ -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; }; diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem_impl.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem_impl.h index b10c6949e6..57ac7d6896 100644 --- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem_impl.h +++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem_impl.h @@ -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; } diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod.h index 61104a143c..48f3301320 100644 --- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod.h +++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod.h @@ -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" diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod1D_impl.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod1D_impl.h index 874945eafd..105f1ceda9 100644 --- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod1D_impl.h +++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod1D_impl.h @@ -55,16 +55,16 @@ template< typename Real, void FastSweepingMethod< Meshes::Grid< 1, Real, Device, Index >, Anisotropy >:: solve( const MeshPointer& mesh, - const AnisotropyType& anisotropy, - MeshFunctionType& u ) + const AnisotropyPointer& anisotropy, + MeshFunctionPointer& u ) { - MeshFunctionType aux; - InterfaceMapType interfaceMap; - aux.setMesh( mesh ); - interfaceMap.setMesh( mesh ); + MeshFunctionPointer aux; + InterfaceMapPointer interfaceMap; + aux->setMesh( mesh ); + interfaceMap->setMesh( mesh ); std::cout << "Initiating the interface cells ..." << std::endl; BaseType::initInterface( u, aux, interfaceMap ); - aux.save( "aux-ini.tnl" ); + aux->save( "aux-ini.tnl" ); } diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod2D_impl.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod2D_impl.h index 48640e3fc5..e36224d075 100644 --- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod2D_impl.h +++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod2D_impl.h @@ -55,22 +55,24 @@ template< typename Real, void FastSweepingMethod< Meshes::Grid< 2, Real, Device, Index >, Anisotropy >:: solve( const MeshPointer& mesh, - const AnisotropyType& anisotropy, - MeshFunctionType& u ) + const AnisotropyPointer& anisotropy, + MeshFunctionPointer& u ) { - MeshFunctionType aux; - InterfaceMapType interfaceMap; - aux.setMesh( mesh ); - interfaceMap.setMesh( mesh ); + MeshFunctionPointer auxPtr; + InterfaceMapPointer interfaceMapPtr; + auxPtr->setMesh( mesh ); + interfaceMapPtr->setMesh( mesh ); std::cout << "Initiating the interface cells ..." << std::endl; - BaseType::initInterface( u, aux, interfaceMap ); + BaseType::initInterface( u, auxPtr, interfaceMapPtr ); cudaDeviceSynchronize(); - aux.save( "aux-ini.tnl" ); + auxPtr->save( "aux-ini.tnl" ); typename MeshType::Cell cell( *mesh ); IndexType iteration( 0 ); + InterfaceMapType interfaceMap = *interfaceMapPtr; + MeshFunctionType aux = *auxPtr; while( iteration < this->maxIterations ) { if( std::is_same< DeviceType, Devices::Host >::value ) diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod3D_impl.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod3D_impl.h index 9844c4fcaf..1fd8c9f2a0 100644 --- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod3D_impl.h +++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod3D_impl.h @@ -55,20 +55,22 @@ template< typename Real, void FastSweepingMethod< Meshes::Grid< 3, Real, Device, Index >, Anisotropy >:: solve( const MeshPointer& mesh, - const AnisotropyType& anisotropy, - MeshFunctionType& u ) + const AnisotropyPointer& anisotropy, + MeshFunctionPointer& u ) { - MeshFunctionType aux; - InterfaceMapType interfaceMap; - aux.setMesh( mesh ); - interfaceMap.setMesh( mesh ); + MeshFunctionPointer auxPtr; + InterfaceMapPointer interfaceMapPtr; + auxPtr->setMesh( mesh ); + interfaceMapPtr->setMesh( mesh ); std::cout << "Initiating the interface cells ..." << std::endl; - BaseType::initInterface( u, aux, interfaceMap ); - aux.save( "aux-ini.tnl" ); + BaseType::initInterface( u, auxPtr, interfaceMapPtr ); + auxPtr->save( "aux-ini.tnl" ); typename MeshType::Cell cell( *mesh ); IndexType iteration( 0 ); + MeshFunctionType aux = *auxPtr; + InterfaceMapType interfaceMap = * interfaceMapPtr; while( iteration < this->maxIterations ) { for( cell.getCoordinates().z() = 0; -- GitLab