Loading src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnl-direct-eikonal-solver.h +1 −1 Original line number Diff line number Diff line Loading @@ -36,7 +36,7 @@ class DirectEikonalSolverConfig { config.addDelimiter( "Direct eikonal equation solver settings:" ); config.addRequiredEntry< String >( "input-file", "Input file." ); config.addEntry< String >( "distributed-grid-io-type", "Choose Distributed Grid IO Type", "LocalCopy"); config.addEntry< String >( "distributed-grid-io-type", "Choose Distributed Grid IO Type", "MpiIO"); config.addEntryEnum< String >( "LocalCopy" ); config.addEntryEnum< String >( "MpiIO" ); }; Loading src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h +1 −1 Original line number Diff line number Diff line Loading @@ -152,7 +152,7 @@ __global__ void CudaUpdateCellCaller( tnlDirectEikonalMethodsBase< Meshes::Grid< template< typename Real, typename Device, typename Index > __global__ void DeepCopy( const Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index > >& aux, Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index > >& helpFunc ); Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index > >& helpFunc, int copy, int k ); template < typename Index > __global__ void CudaParallelReduc( TNL::Containers::ArrayView< int, Devices::Cuda, Index > BlockIterDevice, Loading src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h +34 −69 Original line number Diff line number Diff line Loading @@ -771,8 +771,8 @@ initInterface( const MeshFunctionPointer& _input, { cell.refresh(); output[ cell.getIndex() ] = input( cell ) >= 0 ? 10://std::numeric_limits< RealType >::max() : -10;//- std::numeric_limits< RealType >::max(); input( cell ) >= 0 ? std::numeric_limits< RealType >::max() : - std::numeric_limits< RealType >::max(); interfaceMap[ cell.getIndex() ] = false; } Loading Loading @@ -908,8 +908,8 @@ updateCell( MeshFunctionType& u, b = TNL::argAbsMin( u[ neighborEntities.template getEntityIndex< 0, -1 >() ], u[ neighborEntities.template getEntityIndex< 0, 1 >() ] ); } if( fabs( a ) == 10&&//std::numeric_limits< RealType >::max() && fabs( b ) == 10)//std::numeric_limits< RealType >::max() ) if( fabs( a ) == std::numeric_limits< RealType >::max() && fabs( b ) == std::numeric_limits< RealType >::max() ) return false; RealType pom[6] = { a, b, std::numeric_limits< RealType >::max(), (RealType)hx, (RealType)hy, 0.0 }; Loading Loading @@ -1022,77 +1022,42 @@ initInterface( const MeshFunctionPointer& _input, const IndexType n = neighbors.template getEntityIndex< 0, 1, 0 >(); const IndexType t = neighbors.template getEntityIndex< 0, 0, 1 >(); if( c * input[ n ] <= 0 ) { if( c >= 0 ) { pom = ( hy * c )/( c - input[ n ]); if( output[ cell.getIndex() ] > pom ) output[ cell.getIndex() ] = pom; if ( output[ n ] < pom - hy) output[ n ] = pom - hy; // ( hy * c )/( c - input[ n ]) - hy; }else if( c * input[ n ] <= 0 ) { pom = - ( hy * c )/( c - input[ n ]); if( output[ cell.getIndex() ] < pom ) pom = TNL::sign( c )*( hy * c )/( c - input[ n ]); if( TNL::abs( output[ cell.getIndex() ] ) > TNL::abs( pom ) ) output[ cell.getIndex() ] = pom; if( output[ n ] > hy + pom ) output[ n ] = hy + pom; //hy - ( hy * c )/( c - input[ n ]); pom = pom - TNL::sign( c )*hy; if( TNL::abs( output[ n ] ) > TNL::abs( pom ) ) output[ n ] = pom; //( hy * c )/( c - input[ n ]) - hy; } interfaceMap[ cell.getIndex() ] = true; interfaceMap[ n ] = true; } if( c * input[ e ] <= 0 ) { if( c >= 0 ) { pom = ( hx * c )/( c - input[ e ]); if( output[ cell.getIndex() ] > pom ) output[ cell.getIndex() ] = pom; pom = pom - hx; //output[ e ] = (hx * c)/( c - input[ e ]) - hx; if( output[ e ] < pom ) output[ e ] = pom; }else if( c * input[ e ] <= 0 ) { pom = - (hx * c)/( c - input[ e ]); if( output[ cell.getIndex() ] < pom ) pom = TNL::sign( c )*( hx * c )/( c - input[ e ]); if( TNL::abs( output[ cell.getIndex() ] ) > TNL::abs( pom ) ) output[ cell.getIndex() ] = pom; pom = pom - TNL::sign( c )*hx; if( TNL::abs( output[ e ] ) > TNL::abs( pom ) ) output[ e ] = pom; //( hy * c )/( c - input[ n ]) - hy; pom = pom + hx; //output[ e ] = hx - (hx * c)/( c - input[ e ]); if( output[ e ] > pom ) output[ e ] = pom; } interfaceMap[ cell.getIndex() ] = true; interfaceMap[ e ] = true; } if( c * input[ t ] <= 0 ) { if( c >= 0 ) { pom = ( hz * c )/( c - input[ t ]); if( output[ cell.getIndex() ] > pom ) output[ cell.getIndex() ] = pom; pom = pom - hz; //output[ e ] = (hx * c)/( c - input[ e ]) - hx; if( output[ t ] < pom ) output[ t ] = pom; }else if( c * input[ t ] <= 0 ) { pom = - (hz * c)/( c - input[ t ]); if( output[ cell.getIndex() ] < pom ) pom = TNL::sign( c )*( hz * c )/( c - input[ t ]); if( TNL::abs( output[ cell.getIndex() ] ) > TNL::abs( pom ) ) output[ cell.getIndex() ] = pom; pom = pom - TNL::sign( c )*hz; if( TNL::abs( output[ t ] ) > TNL::abs( pom ) ) output[ t ] = pom; //( hy * c )/( c - input[ n ]) - hy; pom = pom + hz; //output[ e ] = hx - (hx * c)/( c - input[ e ]); if( output[ t ] > pom ) output[ t ] = pom; } interfaceMap[ cell.getIndex() ] = true; interfaceMap[ t ] = true; } Loading Loading @@ -1240,8 +1205,8 @@ updateCell( volatile Real *sArray, int thri, int thrj, const Real hx, const Real a = TNL::argAbsMin( sArray[ thrj * sizeSArray + thri+1 ], sArray[ thrj * sizeSArray + thri-1 ] ); if( fabs( a ) == 10&&//std::numeric_limits< RealType >::max() && fabs( b ) == 10)//std::numeric_limits< RealType >::max() ) if( fabs( a ) == std::numeric_limits< RealType >::max() && fabs( b ) == std::numeric_limits< RealType >::max() ) return false; RealType pom[6] = { a, b, std::numeric_limits< RealType >::max(), (RealType)hx, (RealType)hy, 0.0 }; Loading Loading @@ -1422,8 +1387,8 @@ __global__ void CudaInitCaller( const Functions::MeshFunction< Meshes::Grid< 2, output[ cind ] = input( cell ) >= 0 ? 10://std::numeric_limits< Real >::max() : - 10;//- std::numeric_limits< Real >::max(); input( cell ) >= 0 ? std::numeric_limits< Real >::max() : - std::numeric_limits< Real >::max(); interfaceMap[ cind ] = false; if( i < mesh.getDimensions().x() - vUpper[0] && j < mesh.getDimensions().y() - vUpper[1] && i>vLower[0] && j> vLower[0] ) Loading src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem_impl.h +1 −0 Original line number Diff line number Diff line Loading @@ -191,6 +191,7 @@ bool tnlDirectEikonalProblem< Mesh, Communicator, Anisotropy, Real, Index >:: solve( DofVectorPointer& dofs ) { std::cout << "We are in solve()." << std::endl; FastSweepingMethod< MeshType, Communicator,AnisotropyType > fsm; fsm.solve( this->getMesh(), u, anisotropy, initialData ); Loading src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod.h +2 −0 Original line number Diff line number Diff line Loading @@ -86,6 +86,8 @@ class FastSweepingMethod< Meshes::Grid< 2, Real, Device, Index >, Communicator, typedef Index IndexType; typedef Anisotropy AnisotropyType; typedef tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > > BaseType; typedef Communicator CommunicatorType; using MeshPointer = Pointers::SharedPointer< MeshType >; using AnisotropyPointer = Pointers::SharedPointer< AnisotropyType, DeviceType >; using MPI = Communicators::MpiCommunicator; Loading Loading
src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnl-direct-eikonal-solver.h +1 −1 Original line number Diff line number Diff line Loading @@ -36,7 +36,7 @@ class DirectEikonalSolverConfig { config.addDelimiter( "Direct eikonal equation solver settings:" ); config.addRequiredEntry< String >( "input-file", "Input file." ); config.addEntry< String >( "distributed-grid-io-type", "Choose Distributed Grid IO Type", "LocalCopy"); config.addEntry< String >( "distributed-grid-io-type", "Choose Distributed Grid IO Type", "MpiIO"); config.addEntryEnum< String >( "LocalCopy" ); config.addEntryEnum< String >( "MpiIO" ); }; Loading
src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h +1 −1 Original line number Diff line number Diff line Loading @@ -152,7 +152,7 @@ __global__ void CudaUpdateCellCaller( tnlDirectEikonalMethodsBase< Meshes::Grid< template< typename Real, typename Device, typename Index > __global__ void DeepCopy( const Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index > >& aux, Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index > >& helpFunc ); Functions::MeshFunction< Meshes::Grid< 2, Real, Device, Index > >& helpFunc, int copy, int k ); template < typename Index > __global__ void CudaParallelReduc( TNL::Containers::ArrayView< int, Devices::Cuda, Index > BlockIterDevice, Loading
src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase_impl.h +34 −69 Original line number Diff line number Diff line Loading @@ -771,8 +771,8 @@ initInterface( const MeshFunctionPointer& _input, { cell.refresh(); output[ cell.getIndex() ] = input( cell ) >= 0 ? 10://std::numeric_limits< RealType >::max() : -10;//- std::numeric_limits< RealType >::max(); input( cell ) >= 0 ? std::numeric_limits< RealType >::max() : - std::numeric_limits< RealType >::max(); interfaceMap[ cell.getIndex() ] = false; } Loading Loading @@ -908,8 +908,8 @@ updateCell( MeshFunctionType& u, b = TNL::argAbsMin( u[ neighborEntities.template getEntityIndex< 0, -1 >() ], u[ neighborEntities.template getEntityIndex< 0, 1 >() ] ); } if( fabs( a ) == 10&&//std::numeric_limits< RealType >::max() && fabs( b ) == 10)//std::numeric_limits< RealType >::max() ) if( fabs( a ) == std::numeric_limits< RealType >::max() && fabs( b ) == std::numeric_limits< RealType >::max() ) return false; RealType pom[6] = { a, b, std::numeric_limits< RealType >::max(), (RealType)hx, (RealType)hy, 0.0 }; Loading Loading @@ -1022,77 +1022,42 @@ initInterface( const MeshFunctionPointer& _input, const IndexType n = neighbors.template getEntityIndex< 0, 1, 0 >(); const IndexType t = neighbors.template getEntityIndex< 0, 0, 1 >(); if( c * input[ n ] <= 0 ) { if( c >= 0 ) { pom = ( hy * c )/( c - input[ n ]); if( output[ cell.getIndex() ] > pom ) output[ cell.getIndex() ] = pom; if ( output[ n ] < pom - hy) output[ n ] = pom - hy; // ( hy * c )/( c - input[ n ]) - hy; }else if( c * input[ n ] <= 0 ) { pom = - ( hy * c )/( c - input[ n ]); if( output[ cell.getIndex() ] < pom ) pom = TNL::sign( c )*( hy * c )/( c - input[ n ]); if( TNL::abs( output[ cell.getIndex() ] ) > TNL::abs( pom ) ) output[ cell.getIndex() ] = pom; if( output[ n ] > hy + pom ) output[ n ] = hy + pom; //hy - ( hy * c )/( c - input[ n ]); pom = pom - TNL::sign( c )*hy; if( TNL::abs( output[ n ] ) > TNL::abs( pom ) ) output[ n ] = pom; //( hy * c )/( c - input[ n ]) - hy; } interfaceMap[ cell.getIndex() ] = true; interfaceMap[ n ] = true; } if( c * input[ e ] <= 0 ) { if( c >= 0 ) { pom = ( hx * c )/( c - input[ e ]); if( output[ cell.getIndex() ] > pom ) output[ cell.getIndex() ] = pom; pom = pom - hx; //output[ e ] = (hx * c)/( c - input[ e ]) - hx; if( output[ e ] < pom ) output[ e ] = pom; }else if( c * input[ e ] <= 0 ) { pom = - (hx * c)/( c - input[ e ]); if( output[ cell.getIndex() ] < pom ) pom = TNL::sign( c )*( hx * c )/( c - input[ e ]); if( TNL::abs( output[ cell.getIndex() ] ) > TNL::abs( pom ) ) output[ cell.getIndex() ] = pom; pom = pom - TNL::sign( c )*hx; if( TNL::abs( output[ e ] ) > TNL::abs( pom ) ) output[ e ] = pom; //( hy * c )/( c - input[ n ]) - hy; pom = pom + hx; //output[ e ] = hx - (hx * c)/( c - input[ e ]); if( output[ e ] > pom ) output[ e ] = pom; } interfaceMap[ cell.getIndex() ] = true; interfaceMap[ e ] = true; } if( c * input[ t ] <= 0 ) { if( c >= 0 ) { pom = ( hz * c )/( c - input[ t ]); if( output[ cell.getIndex() ] > pom ) output[ cell.getIndex() ] = pom; pom = pom - hz; //output[ e ] = (hx * c)/( c - input[ e ]) - hx; if( output[ t ] < pom ) output[ t ] = pom; }else if( c * input[ t ] <= 0 ) { pom = - (hz * c)/( c - input[ t ]); if( output[ cell.getIndex() ] < pom ) pom = TNL::sign( c )*( hz * c )/( c - input[ t ]); if( TNL::abs( output[ cell.getIndex() ] ) > TNL::abs( pom ) ) output[ cell.getIndex() ] = pom; pom = pom - TNL::sign( c )*hz; if( TNL::abs( output[ t ] ) > TNL::abs( pom ) ) output[ t ] = pom; //( hy * c )/( c - input[ n ]) - hy; pom = pom + hz; //output[ e ] = hx - (hx * c)/( c - input[ e ]); if( output[ t ] > pom ) output[ t ] = pom; } interfaceMap[ cell.getIndex() ] = true; interfaceMap[ t ] = true; } Loading Loading @@ -1240,8 +1205,8 @@ updateCell( volatile Real *sArray, int thri, int thrj, const Real hx, const Real a = TNL::argAbsMin( sArray[ thrj * sizeSArray + thri+1 ], sArray[ thrj * sizeSArray + thri-1 ] ); if( fabs( a ) == 10&&//std::numeric_limits< RealType >::max() && fabs( b ) == 10)//std::numeric_limits< RealType >::max() ) if( fabs( a ) == std::numeric_limits< RealType >::max() && fabs( b ) == std::numeric_limits< RealType >::max() ) return false; RealType pom[6] = { a, b, std::numeric_limits< RealType >::max(), (RealType)hx, (RealType)hy, 0.0 }; Loading Loading @@ -1422,8 +1387,8 @@ __global__ void CudaInitCaller( const Functions::MeshFunction< Meshes::Grid< 2, output[ cind ] = input( cell ) >= 0 ? 10://std::numeric_limits< Real >::max() : - 10;//- std::numeric_limits< Real >::max(); input( cell ) >= 0 ? std::numeric_limits< Real >::max() : - std::numeric_limits< Real >::max(); interfaceMap[ cind ] = false; if( i < mesh.getDimensions().x() - vUpper[0] && j < mesh.getDimensions().y() - vUpper[1] && i>vLower[0] && j> vLower[0] ) Loading
src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem_impl.h +1 −0 Original line number Diff line number Diff line Loading @@ -191,6 +191,7 @@ bool tnlDirectEikonalProblem< Mesh, Communicator, Anisotropy, Real, Index >:: solve( DofVectorPointer& dofs ) { std::cout << "We are in solve()." << std::endl; FastSweepingMethod< MeshType, Communicator,AnisotropyType > fsm; fsm.solve( this->getMesh(), u, anisotropy, initialData ); Loading
src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlFastSweepingMethod.h +2 −0 Original line number Diff line number Diff line Loading @@ -86,6 +86,8 @@ class FastSweepingMethod< Meshes::Grid< 2, Real, Device, Index >, Communicator, typedef Index IndexType; typedef Anisotropy AnisotropyType; typedef tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > > BaseType; typedef Communicator CommunicatorType; using MeshPointer = Pointers::SharedPointer< MeshType >; using AnisotropyPointer = Pointers::SharedPointer< AnisotropyType, DeviceType >; using MPI = Communicators::MpiCommunicator; Loading