From d2d59429c4901bd8490020f25a7c2f160c12a78a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matou=C5=A1=20Fencl?= Date: Fri, 30 Nov 2018 09:02:24 +0100 Subject: [PATCH 1/6] Change of int to IndexType and preparations for OpenMPI. --- .../hamilton-jacobi/HamiltonJacobiProblem.h | 2 ++ .../HamiltonJacobiProblem_impl.h | 24 ++++++++++++++----- 2 files changed, 20 insertions(+), 6 deletions(-) diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem.h index a41442000..7f1bd4193 100644 --- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem.h +++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem.h @@ -25,6 +25,8 @@ #include #include +#include + template< typename Mesh, typename DifferentialOperator, typename BoundaryCondition, diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem_impl.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem_impl.h index b529b09fd..9244b1833 100644 --- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem_impl.h +++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem_impl.h @@ -121,13 +121,25 @@ setInitialCondition( const Config::ParameterContainer& parameters, DofVectorType& dofs, MeshDependentDataType& meshDependentData ) { - this->bindDofs( mesh, dofs ); - const String& initialConditionFile = parameters.getParameter< String >( "initial-condition" ); - if( ! this->solution.boundLoad( initialConditionFile ) ) - { - std::cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." <bindDofs( mesh, dofs ); + const String& initialConditionFile = parameters.getParameter< String >( "initial-condition" ); + if(CommunicatorType::isDistributed()) + { + std::cout<<"Nodes Distribution: " << uPointer->getMesh().getDistributedMesh()->printProcessDistr() << std::endl; + if(distributedIOType==Meshes::DistributedMeshes::MpiIO) + Meshes::DistributedMeshes::DistributedGridIO ::load(initialConditionFile, *uPointer ); + if(distributedIOType==Meshes::DistributedMeshes::LocalCopy) + Meshes::DistributedMeshes::DistributedGridIO ::load(initialConditionFile, *uPointer ); + uPointer->template synchronize(); + } + else + { + if( ! this->solution.boundLoad( initialConditionFile ) ) + { + std::cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." < Date: Fri, 30 Nov 2018 09:06:24 +0100 Subject: [PATCH 2/6] Changed int to IndexType --- .../tnlDirectEikonalMethodsBase.h | 3 - .../tnlDirectEikonalMethodsBase_impl.h | 132 ++++++++---------- .../tnlFastSweepingMethod2D_impl.h | 29 +++- 3 files changed, 81 insertions(+), 83 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 f712ce2cc..76d98a2de 100644 --- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h +++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalMethodsBase.h @@ -127,9 +127,6 @@ class tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > > const Real velocity = 1.0 ); }; -template < typename T1, typename T2 > -T1 meet2DCondition( T1 a, T1 b, const T2 ha, const T2 hb, const T1 value, double v = 1); - template < typename T1 > __cuda_callable__ void sortMinims( T1 pom[] ); 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 8f7937541..e142152ef 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 @@ -99,10 +99,10 @@ tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > >:: updateBlocks( InterfaceMapType interfaceMap, MeshFunctionType aux, MeshFunctionType helpFunc, - ArrayContainer BlockIterHost, int numThreadsPerBlock/*, Real **sArray*/ ) + ArrayContainer BlockIterHost, IndexType numThreadsPerBlock/*, Real **sArray*/ ) { #pragma omp parallel for schedule( dynamic ) - for( int i = 0; i < BlockIterHost.getSize(); i++ ) + for( IndexType i = 0; i < BlockIterHost.getSize(); i++ ) { if( BlockIterHost[ i ] ) { @@ -142,14 +142,14 @@ updateBlocks( InterfaceMapType interfaceMap, if( sArray == nullptr ) std::cout << "Error while allocating memory for sArray." << std::endl; - for( int thri = 0; thri < sizeSArray; thri++ ){ - for( int thrj = 0; thrj < sizeSArray; thrj++ ) + for( IndexType thri = 0; thri < sizeSArray; thri++ ){ + for( IndexType thrj = 0; thrj < sizeSArray; thrj++ ) sArray[ thri * sizeSArray + thrj ] = std::numeric_limits< Real >::max(); } //printf("numThreadsPerBlock = %d\n", numThreadsPerBlock); - for( int thrj = 0; thrj < numThreadsPerBlock + 1; thrj++ ) + for( IndexType thrj = 0; thrj < numThreadsPerBlock + 1; thrj++ ) { if( dimX > (blIdx+1) * numThreadsPerBlock && thrj+1 < ykolik ) sArray[ ( thrj+1 )* sizeSArray +xkolik] = aux[ blIdy*numThreadsPerBlock*dimX - dimX + blIdx*numThreadsPerBlock - 1 + (thrj+1)*dimX + xkolik ]; @@ -165,14 +165,14 @@ updateBlocks( InterfaceMapType interfaceMap, sArray[thrj+1] = aux[ blIdy*numThreadsPerBlock*dimX - dimX + blIdx*numThreadsPerBlock - 1 + thrj+1 ]; } - for( int k = 0; k < numThreadsPerBlock; k++ ){ - for( int l = 0; l < numThreadsPerBlock; l++ ) + for( IndexType k = 0; k < numThreadsPerBlock; k++ ){ + for( IndexType l = 0; l < numThreadsPerBlock; l++ ) if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX ) sArray[(k+1) * sizeSArray + l+1] = aux[ blIdy * numThreadsPerBlock * dimX + numThreadsPerBlock * blIdx + k*dimX + l ]; } - for( int k = 0; k < numThreadsPerBlock; k++ ){ - for( int l = 0; l < numThreadsPerBlock; l++ ){ + for( IndexType k = 0; k < numThreadsPerBlock; k++ ){ + for( IndexType l = 0; l < numThreadsPerBlock; l++ ){ if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX ){ //std::cout << "proslo i = " << k * numThreadsPerBlock + l << std::endl; if( ! interfaceMap[ blIdy * numThreadsPerBlock * dimX + numThreadsPerBlock * blIdx + k*dimX + l ] ) @@ -191,8 +191,8 @@ updateBlocks( InterfaceMapType interfaceMap, std::cout << std::endl; }*/ - for( int k = 0; k < numThreadsPerBlock; k++ ) - for( int l = numThreadsPerBlock-1; l >-1; l-- ) { + for( IndexType k = 0; k < numThreadsPerBlock; k++ ) + for( IndexType l = numThreadsPerBlock-1; l >-1; l-- ) { if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX ) { if( ! interfaceMap[ blIdy * numThreadsPerBlock * dimX + numThreadsPerBlock * blIdx + k*dimX + l ] ) @@ -209,8 +209,8 @@ updateBlocks( InterfaceMapType interfaceMap, std::cout << std::endl; }*/ - for( int k = numThreadsPerBlock-1; k > -1; k-- ) - for( int l = 0; l < numThreadsPerBlock; l++ ) { + for( IndexType k = numThreadsPerBlock-1; k > -1; k-- ) + for( IndexType l = 0; l < numThreadsPerBlock; l++ ) { if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX ) { if( ! interfaceMap[ blIdy * numThreadsPerBlock * dimX + numThreadsPerBlock * blIdx + k*dimX + l ] ) @@ -227,8 +227,8 @@ updateBlocks( InterfaceMapType interfaceMap, std::cout << std::endl; }*/ - for( int k = numThreadsPerBlock-1; k > -1; k-- ){ - for( int l = numThreadsPerBlock-1; l >-1; l-- ) { + for( IndexType k = numThreadsPerBlock-1; k > -1; k-- ){ + for( IndexType l = numThreadsPerBlock-1; l >-1; l-- ) { if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX ) { if( ! interfaceMap[ blIdy * numThreadsPerBlock * dimX + numThreadsPerBlock * blIdx + k*dimX + l ] ) @@ -252,8 +252,8 @@ updateBlocks( InterfaceMapType interfaceMap, } - for( int k = 0; k < numThreadsPerBlock; k++ ){ - for( int l = 0; l < numThreadsPerBlock; l++ ) { + for( IndexType k = 0; k < numThreadsPerBlock; k++ ){ + for( IndexType l = 0; l < numThreadsPerBlock; l++ ) { if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX ) helpFunc[ blIdy * numThreadsPerBlock * dimX + numThreadsPerBlock * blIdx + k*dimX + l ] = sArray[ (k + 1)* sizeSArray + l + 1 ]; //std::cout<< sArray[k+1][l+1]; @@ -267,16 +267,16 @@ updateBlocks( InterfaceMapType interfaceMap, template< typename Real, typename Device, typename Index > -template< int sizeSArray > +template< IndexType sizeSArray > void tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > >:: updateBlocks( const InterfaceMapType interfaceMap, const MeshFunctionType aux, MeshFunctionType& helpFunc, - ArrayContainer BlockIterHost, int numThreadsPerBlock/*, Real **sArray*/ ) + ArrayContainer BlockIterHost, IndexType numThreadsPerBlock/*, Real **sArray*/ ) { //#pragma omp parallel for schedule( dynamic ) - for( int i = 0; i < BlockIterHost.getSize(); i++ ) + for( IndexType i = 0; i < BlockIterHost.getSize(); i++ ) { if( BlockIterHost[ i ] ) { @@ -323,15 +323,15 @@ updateBlocks( const InterfaceMapType interfaceMap, if( sArray == nullptr ) std::cout << "Error while allocating memory for sArray." << std::endl; - for( int k = 0; k < sizeSArray; k++ ) - for( int l = 0; l < sizeSArray; l++ ) - for( int m = 0; m < sizeSArray; m++ ){ + for( IndexType k = 0; k < sizeSArray; k++ ) + for( IndexType l = 0; l < sizeSArray; l++ ) + for( IndexType m = 0; m < sizeSArray; m++ ){ sArray[ m * sizeSArray * sizeSArray + k * sizeSArray + l ] = std::numeric_limits< Real >::max(); } - for( int thrk = 0; thrk < numThreadsPerBlock; thrk++ ) - for( int thrj = 0; thrj < numThreadsPerBlock; thrj++ ) + for( IndexType thrk = 0; thrk < numThreadsPerBlock; thrk++ ) + for( IndexType thrj = 0; thrj < numThreadsPerBlock; thrj++ ) { if( blIdx != 0 && thrj+1 < ykolik && thrk+1 < zkolik ) sArray[(thrk+1 )* sizeSArray * sizeSArray + (thrj+1)*sizeSArray + 0] = @@ -358,9 +358,9 @@ updateBlocks( const InterfaceMapType interfaceMap, aux[ (blIdz+1)*numThreadsPerBlock * dimX * dimY + blIdy * numThreadsPerBlock*dimX + blIdx*numThreadsPerBlock + thrj * dimX + thrk ]; } - for( int m = 0; m < numThreadsPerBlock; m++ ){ - for( int k = 0; k < numThreadsPerBlock; k++ ){ - for( int l = 0; l < numThreadsPerBlock; l++ ){ + for( IndexType m = 0; m < numThreadsPerBlock; m++ ){ + for( IndexType k = 0; k < numThreadsPerBlock; k++ ){ + for( IndexType l = 0; l < numThreadsPerBlock; l++ ){ if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX && blIdz * numThreadsPerBlock + m < dimZ ) sArray[(m+1) * sizeSArray * sizeSArray + (k+1) *sizeSArray + l+1] = aux[ blIdz * numThreadsPerBlock * dimX * dimY + blIdy * numThreadsPerBlock * dimX + blIdx*numThreadsPerBlock + m*dimX*dimY + k*dimX + l ]; @@ -379,9 +379,9 @@ updateBlocks( const InterfaceMapType interfaceMap, s = "helpFunc-"+ std::to_string(numWhile) + ".tnl"; helpFunc.save( s );*/ - for( int m = 0; m < numThreadsPerBlock; m++ ){ - for( int k = 0; k < numThreadsPerBlock; k++ ){ - for( int l = 0; l < numThreadsPerBlock; l++ ){ + for( IndexType m = 0; m < numThreadsPerBlock; m++ ){ + for( IndexType k = 0; k < numThreadsPerBlock; k++ ){ + for( IndexType l = 0; l < numThreadsPerBlock; l++ ){ if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX && blIdz * numThreadsPerBlock + m < dimZ ){ //std::cout << "proslo i = " << k * numThreadsPerBlock + l << std::endl; if( ! interfaceMap[ blIdz * numThreadsPerBlock * dimX * dimY + blIdy * numThreadsPerBlock * dimX + blIdx*numThreadsPerBlock + m*dimX*dimY + k*dimX + l ] ) @@ -404,9 +404,9 @@ updateBlocks( const InterfaceMapType interfaceMap, s = "helpFunc-"+ std::to_string(numWhile) + ".tnl"; helpFunc.save( s );*/ - for( int m = numThreadsPerBlock-1; m >-1; m-- ){ - for( int k = 0; k < numThreadsPerBlock; k++ ){ - for( int l = 0; l -1; m-- ){ + for( IndexType k = 0; k < numThreadsPerBlock; k++ ){ + for( IndexType l = 0; l -1; l-- ){ + for( IndexType m = 0; m < numThreadsPerBlock; m++ ){ + for( IndexType k = 0; k < numThreadsPerBlock; k++ ){ + for( IndexType l = numThreadsPerBlock-1; l >-1; l-- ){ if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX && blIdz * numThreadsPerBlock + m < dimZ ) { if( ! interfaceMap[ blIdz * numThreadsPerBlock * dimX * dimY + blIdy * numThreadsPerBlock * dimX + blIdx*numThreadsPerBlock + m*dimX*dimY + k*dimX + l ] ) @@ -450,9 +450,9 @@ updateBlocks( const InterfaceMapType interfaceMap, s = "helpFunc-"+ std::to_string(numWhile) + ".tnl"; helpFunc.save( s ); */ - for( int m = numThreadsPerBlock-1; m >-1; m-- ){ - for( int k = 0; k < numThreadsPerBlock; k++ ){ - for( int l = numThreadsPerBlock-1; l >-1; l-- ){ + for( IndexType m = numThreadsPerBlock-1; m >-1; m-- ){ + for( IndexType k = 0; k < numThreadsPerBlock; k++ ){ + for( IndexType l = numThreadsPerBlock-1; l >-1; l-- ){ if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX && blIdz * numThreadsPerBlock + m < dimZ ) { if( ! interfaceMap[ blIdz * numThreadsPerBlock * dimX * dimY + blIdy * numThreadsPerBlock * dimX + blIdx*numThreadsPerBlock + m*dimX*dimY + k*dimX + l ] ) @@ -473,9 +473,9 @@ updateBlocks( const InterfaceMapType interfaceMap, s = "helpFunc-"+ std::to_string(numWhile) + ".tnl"; helpFunc.save( s );*/ - for( int m = 0; m < numThreadsPerBlock; m++ ){ - for( int k = numThreadsPerBlock-1; k > -1; k-- ){ - for( int l = 0; l -1; k-- ){ + for( IndexType l = 0; l -1; m-- ){ - for( int k = numThreadsPerBlock-1; k > -1; k-- ){ - for( int l = 0; l -1; m-- ){ + for( IndexType k = numThreadsPerBlock-1; k > -1; k-- ){ + for( IndexType l = 0; l -1; k-- ){ - for( int l = numThreadsPerBlock-1; l >-1; l-- ){ + for( IndexType m = 0; m < numThreadsPerBlock; m++ ){ + for( IndexType k = numThreadsPerBlock-1; k > -1; k-- ){ + for( IndexType l = numThreadsPerBlock-1; l >-1; l-- ){ if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX && blIdz * numThreadsPerBlock + m < dimZ ) { if( ! interfaceMap[ blIdz * numThreadsPerBlock * dimX * dimY + blIdy * numThreadsPerBlock * dimX + blIdx*numThreadsPerBlock + m*dimX*dimY + k*dimX + l ] ) @@ -543,9 +543,9 @@ updateBlocks( const InterfaceMapType interfaceMap, helpFunc.save( s );*/ - for( int m = numThreadsPerBlock-1; m >-1; m-- ){ - for( int k = numThreadsPerBlock-1; k > -1; k-- ){ - for( int l = numThreadsPerBlock-1; l >-1; l-- ){ + for( IndexType m = numThreadsPerBlock-1; m >-1; m-- ){ + for( IndexType k = numThreadsPerBlock-1; k > -1; k-- ){ + for( IndexType l = numThreadsPerBlock-1; l >-1; l-- ){ if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX && blIdz * numThreadsPerBlock + m < dimZ ) { if( ! interfaceMap[ blIdz * numThreadsPerBlock * dimX * dimY + blIdy * numThreadsPerBlock * dimX + blIdx*numThreadsPerBlock + m*dimX*dimY + k*dimX + l ] ) @@ -571,9 +571,9 @@ updateBlocks( const InterfaceMapType interfaceMap, } - for( int k = 0; k < numThreadsPerBlock; k++ ){ - for( int l = 0; l < numThreadsPerBlock; l++ ) { - for( int m = 0; m < numThreadsPerBlock; m++ ){ + for( IndexType k = 0; k < numThreadsPerBlock; k++ ){ + for( IndexType l = 0; l < numThreadsPerBlock; l++ ) { + for( IndexType m = 0; m < numThreadsPerBlock; m++ ){ if( blIdy * numThreadsPerBlock + k < dimY && blIdx * numThreadsPerBlock + l < dimX && blIdz * numThreadsPerBlock + m < dimZ ){ helpFunc[ blIdz * numThreadsPerBlock * dimX * dimY + blIdy * numThreadsPerBlock * dimX + blIdx*numThreadsPerBlock + m*dimX*dimY + k*dimX + l ] = sArray[ (m+1) * sizeSArray * sizeSArray + (k+1) *sizeSArray + l+1 ]; @@ -1224,26 +1224,6 @@ updateCell( MeshFunctionType& u, } } -template < typename T1, typename T2 > -T1 meet2DCondition( T1 a, T1 b, const T2 ha, const T2 hb, const T1 value, double v) -{ - T1 tmp; - if( fabs( a ) != std::numeric_limits< T1 >::max && - fabs( b ) != std::numeric_limits< T1 >::max && - fabs( a - b ) < ha/v )//TNL::sqrt( (ha * ha + hb * hb)/2 )/v ) - { - tmp = ( ha * ha * b + hb * hb * a + - TNL::sign( value ) * ha * hb * TNL::sqrt( ( ha * ha + hb * hb )/( v * v ) - - ( a - b ) * ( a - b ) ) )/( ha * ha + hb * hb ); - } - else - { - tmp = std::numeric_limits< T1 >::max; - } - - return tmp; -} - template < typename T1 > __cuda_callable__ void sortMinims( T1 pom[] ) { @@ -1274,7 +1254,7 @@ __cuda_callable__ void sortMinims( T1 pom[] ) tmp[3] = pom[5]; tmp[4] = pom[4]; tmp[5] = pom[3]; } - for( int i = 0; i < 6; i++ ) + for( unsigned int i = 0; i < 6; i++ ) { pom[ i ] = tmp[ i ]; } 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 d5ce1efe1..27d266341 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 @@ -83,13 +83,34 @@ solve( const MeshPointer& mesh, MeshFunctionType aux = *auxPtr; -//#ifdef HAVE_MPI - bool a = Communicators::MpiCommunicator::IsInitialized(); +#ifdef HAVE_MPI + int i = Communicators::MpiCommunicator::GetRank( Communicators::MpiCommunicator::AllGroup ); + printf( "Hello world from rank: %d ", i ); + //Communicators::MpiCommunicator::Request r = Communicators::MpiCommunicator::ISend( auxPtr, 0, 0, Communicators::MpiCommunicator::AllGroup ); + /*if( i == 1 ) + for( int k = 0; k < 16*16; k++ ) + aux[ k ] = 10; + printf( "1: mesh x: %d\n", mesh->getDimensions().x() ); + printf( "1: mesh y: %d\n", mesh->getDimensions().y() ); + aux.save("aux_proc1.tnl");*/ + if( i == 0 ) + printf( "0: mesh x: %d\n", mesh->getDimensions().x() ); + printf( "0: mesh y: %d\n", mesh->getDimensions().y() ); + aux.save("aux_proc0.tnl"); + for( int k = 0; k < 16*16; k++ ) + aux[ k ] = 10; + for( int k = 0; k < 16; k++ ){ + for( int l = 0; l < 16; l++ ) + printf("%f.2\t",aux[ k * 16 + l ] ); + printf("\n"); + } + + /*bool a = Communicators::MpiCommunicator::IsInitialized(); if( a ) printf("Je Init\n"); else - printf("Neni Init\n"); -//#endif + printf("Neni Init\n");*/ +#endif while( iteration < this->maxIterations ) { -- GitLab From b387127c77aed9f619fc2a5439bd31d21772b0ae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matou=C5=A1=20Fencl?= Date: Fri, 30 Nov 2018 15:16:19 +0100 Subject: [PATCH 3/6] MPI ready in tnlDirectEikonal* --- .../hamilton-jacobi/HamiltonJacobiProblem.h | 2 -- .../HamiltonJacobiProblem_impl.h | 20 +++------------ .../tnlDirectEikonalMethodsBase_impl.h | 6 ++--- .../hamilton-jacobi/tnlDirectEikonalProblem.h | 8 +++++- .../tnlDirectEikonalProblem_impl.h | 25 +++++++++++++++---- .../tnlFastSweepingMethod2D_impl.h | 20 +++++++-------- 6 files changed, 44 insertions(+), 37 deletions(-) diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem.h index 7f1bd4193..a41442000 100644 --- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem.h +++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem.h @@ -25,8 +25,6 @@ #include #include -#include - template< typename Mesh, typename DifferentialOperator, typename BoundaryCondition, diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem_impl.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem_impl.h index 9244b1833..3cc638849 100644 --- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem_impl.h +++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem_impl.h @@ -123,24 +123,12 @@ setInitialCondition( const Config::ParameterContainer& parameters, { this->bindDofs( mesh, dofs ); const String& initialConditionFile = parameters.getParameter< String >( "initial-condition" ); - if(CommunicatorType::isDistributed()) + if( ! this->solution.boundLoad( initialConditionFile ) ) { - std::cout<<"Nodes Distribution: " << uPointer->getMesh().getDistributedMesh()->printProcessDistr() << std::endl; - if(distributedIOType==Meshes::DistributedMeshes::MpiIO) - Meshes::DistributedMeshes::DistributedGridIO ::load(initialConditionFile, *uPointer ); - if(distributedIOType==Meshes::DistributedMeshes::LocalCopy) - Meshes::DistributedMeshes::DistributedGridIO ::load(initialConditionFile, *uPointer ); - uPointer->template synchronize(); + std::cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." <solution.boundLoad( initialConditionFile ) ) - { - std::cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." < >:: updateBlocks( InterfaceMapType interfaceMap, MeshFunctionType aux, MeshFunctionType helpFunc, - ArrayContainer BlockIterHost, IndexType numThreadsPerBlock/*, Real **sArray*/ ) + ArrayContainer BlockIterHost, int numThreadsPerBlock/*, Real **sArray*/ ) { #pragma omp parallel for schedule( dynamic ) for( IndexType i = 0; i < BlockIterHost.getSize(); i++ ) @@ -267,13 +267,13 @@ updateBlocks( InterfaceMapType interfaceMap, template< typename Real, typename Device, typename Index > -template< IndexType sizeSArray > +template< int sizeSArray > void tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > >:: updateBlocks( const InterfaceMapType interfaceMap, const MeshFunctionType aux, MeshFunctionType& helpFunc, - ArrayContainer BlockIterHost, IndexType numThreadsPerBlock/*, Real **sArray*/ ) + ArrayContainer BlockIterHost, int numThreadsPerBlock/*, Real **sArray*/ ) { //#pragma omp parallel for schedule( dynamic ) for( IndexType i = 0; i < BlockIterHost.getSize(); i++ ) 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 0dde93802..61465bee9 100644 --- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem.h +++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/tnlDirectEikonalProblem.h @@ -18,6 +18,8 @@ #include #include "tnlFastSweepingMethod.h" +#include + template< typename Mesh, typename Communicator, typename Anisotropy, @@ -46,6 +48,8 @@ class tnlDirectEikonalProblem using MeshPointer = Pointers::SharedPointer< MeshType >; using DofVectorPointer = Pointers::SharedPointer< DofVectorType >; + typedef Communicator CommunicatorType; + static constexpr bool isTimeDependent() { return false; }; static String getType(); @@ -79,6 +83,8 @@ class tnlDirectEikonalProblem AnisotropyPointer anisotropy; -}; + Meshes::DistributedMeshes::DistrGridIOTypes distributedIOType; + + }; #include "tnlDirectEikonalProblem_impl.h" 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 3d5fa15dc..7c65ef94d 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 @@ -113,11 +113,26 @@ tnlDirectEikonalProblem< Mesh, Communicator, Anisotropy, Real, Index >:: setInitialCondition( const Config::ParameterContainer& parameters, DofVectorPointer& dofs ) { - String inputFile = parameters.getParameter< String >( "input-file" ); - this->initialData->setMesh( this->getMesh() ); - if( !this->initialData->boundLoad( inputFile ) ) - return false; - return true; + this->bindDofs( dofs ); + String inputFile = parameters.getParameter< String >( "input-file" ); + this->initialData->setMesh( this->getMesh() ); + std::cout<<"setInitialCondition" <getMesh().getDistributedMesh()->printProcessDistr() << std::endl; + if(distributedIOType==Meshes::DistributedMeshes::MpiIO) + Meshes::DistributedMeshes::DistributedGridIO ::load(inputFile, *u ); + if(distributedIOType==Meshes::DistributedMeshes::LocalCopy) + Meshes::DistributedMeshes::DistributedGridIO ::load(inputFile, *u ); + u->template synchronize(); + } + else + { + if( !this->initialData->boundLoad( inputFile ) ) + std::cerr << "I am not able to load the initial condition from the file " << inputFile << "." << std::endl; + return false; + } + return true; } 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 27d266341..001f4c69a 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 @@ -85,25 +85,25 @@ solve( const MeshPointer& mesh, #ifdef HAVE_MPI int i = Communicators::MpiCommunicator::GetRank( Communicators::MpiCommunicator::AllGroup ); - printf( "Hello world from rank: %d ", i ); + //printf( "Hello world from rank: %d ", i ); //Communicators::MpiCommunicator::Request r = Communicators::MpiCommunicator::ISend( auxPtr, 0, 0, Communicators::MpiCommunicator::AllGroup ); - /*if( i == 1 ) - for( int k = 0; k < 16*16; k++ ) - aux[ k ] = 10; + if( i == 1 ) + /*for( int k = 0; k < 16*16; k++ ) + aux[ k ] = 10;*/ printf( "1: mesh x: %d\n", mesh->getDimensions().x() ); printf( "1: mesh y: %d\n", mesh->getDimensions().y() ); - aux.save("aux_proc1.tnl");*/ + //aux.save("aux_proc1.tnl"); if( i == 0 ) printf( "0: mesh x: %d\n", mesh->getDimensions().x() ); printf( "0: mesh y: %d\n", mesh->getDimensions().y() ); - aux.save("aux_proc0.tnl"); - for( int k = 0; k < 16*16; k++ ) + //aux.save("aux_proc0.tnl"); + /*for( int k = 0; k < mesh->getDimensions().x()*mesh->getDimensions().y(); k++ ) aux[ k ] = 10; - for( int k = 0; k < 16; k++ ){ - for( int l = 0; l < 16; l++ ) + for( int k = 0; k < mesh->getDimensions().x(); k++ ){ + for( int l = 0; l < mesh->getDimensions().y(); l++ ) printf("%f.2\t",aux[ k * 16 + l ] ); printf("\n"); - } + }*/ /*bool a = Communicators::MpiCommunicator::IsInitialized(); if( a ) -- GitLab From 26e178928ace7c0fa371fc5a97987f1c5a517e39 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matou=C5=A1=20Fencl?= Date: Fri, 30 Nov 2018 09:02:24 +0100 Subject: [PATCH 4/6] Change of int to IndexType and preparations for OpenMPI. --- .../hamilton-jacobi/HamiltonJacobiProblem.h | 2 ++ .../HamiltonJacobiProblem_impl.h | 20 +++++++++++++++---- 2 files changed, 18 insertions(+), 4 deletions(-) diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem.h index a41442000..7f1bd4193 100644 --- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem.h +++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem.h @@ -25,6 +25,8 @@ #include #include +#include + template< typename Mesh, typename DifferentialOperator, typename BoundaryCondition, diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem_impl.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem_impl.h index 3cc638849..9244b1833 100644 --- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem_impl.h +++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem_impl.h @@ -123,12 +123,24 @@ setInitialCondition( const Config::ParameterContainer& parameters, { this->bindDofs( mesh, dofs ); const String& initialConditionFile = parameters.getParameter< String >( "initial-condition" ); - if( ! this->solution.boundLoad( initialConditionFile ) ) + if(CommunicatorType::isDistributed()) { - std::cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." <getMesh().getDistributedMesh()->printProcessDistr() << std::endl; + if(distributedIOType==Meshes::DistributedMeshes::MpiIO) + Meshes::DistributedMeshes::DistributedGridIO ::load(initialConditionFile, *uPointer ); + if(distributedIOType==Meshes::DistributedMeshes::LocalCopy) + Meshes::DistributedMeshes::DistributedGridIO ::load(initialConditionFile, *uPointer ); + uPointer->template synchronize(); } - return true; + else + { + if( ! this->solution.boundLoad( initialConditionFile ) ) + { + std::cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." < Date: Fri, 30 Nov 2018 09:06:24 +0100 Subject: [PATCH 5/6] Changed int to IndexType --- .../tnlDirectEikonalMethodsBase_impl.h | 6 +++--- .../tnlFastSweepingMethod2D_impl.h | 20 +++++++++---------- 2 files changed, 13 insertions(+), 13 deletions(-) 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 786c6a488..e142152ef 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 @@ -99,7 +99,7 @@ tnlDirectEikonalMethodsBase< Meshes::Grid< 2, Real, Device, Index > >:: updateBlocks( InterfaceMapType interfaceMap, MeshFunctionType aux, MeshFunctionType helpFunc, - ArrayContainer BlockIterHost, int numThreadsPerBlock/*, Real **sArray*/ ) + ArrayContainer BlockIterHost, IndexType numThreadsPerBlock/*, Real **sArray*/ ) { #pragma omp parallel for schedule( dynamic ) for( IndexType i = 0; i < BlockIterHost.getSize(); i++ ) @@ -267,13 +267,13 @@ updateBlocks( InterfaceMapType interfaceMap, template< typename Real, typename Device, typename Index > -template< int sizeSArray > +template< IndexType sizeSArray > void tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > >:: updateBlocks( const InterfaceMapType interfaceMap, const MeshFunctionType aux, MeshFunctionType& helpFunc, - ArrayContainer BlockIterHost, int numThreadsPerBlock/*, Real **sArray*/ ) + ArrayContainer BlockIterHost, IndexType numThreadsPerBlock/*, Real **sArray*/ ) { //#pragma omp parallel for schedule( dynamic ) for( IndexType i = 0; i < BlockIterHost.getSize(); i++ ) 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 001f4c69a..27d266341 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 @@ -85,25 +85,25 @@ solve( const MeshPointer& mesh, #ifdef HAVE_MPI int i = Communicators::MpiCommunicator::GetRank( Communicators::MpiCommunicator::AllGroup ); - //printf( "Hello world from rank: %d ", i ); + printf( "Hello world from rank: %d ", i ); //Communicators::MpiCommunicator::Request r = Communicators::MpiCommunicator::ISend( auxPtr, 0, 0, Communicators::MpiCommunicator::AllGroup ); - if( i == 1 ) - /*for( int k = 0; k < 16*16; k++ ) - aux[ k ] = 10;*/ + /*if( i == 1 ) + for( int k = 0; k < 16*16; k++ ) + aux[ k ] = 10; printf( "1: mesh x: %d\n", mesh->getDimensions().x() ); printf( "1: mesh y: %d\n", mesh->getDimensions().y() ); - //aux.save("aux_proc1.tnl"); + aux.save("aux_proc1.tnl");*/ if( i == 0 ) printf( "0: mesh x: %d\n", mesh->getDimensions().x() ); printf( "0: mesh y: %d\n", mesh->getDimensions().y() ); - //aux.save("aux_proc0.tnl"); - /*for( int k = 0; k < mesh->getDimensions().x()*mesh->getDimensions().y(); k++ ) + aux.save("aux_proc0.tnl"); + for( int k = 0; k < 16*16; k++ ) aux[ k ] = 10; - for( int k = 0; k < mesh->getDimensions().x(); k++ ){ - for( int l = 0; l < mesh->getDimensions().y(); l++ ) + for( int k = 0; k < 16; k++ ){ + for( int l = 0; l < 16; l++ ) printf("%f.2\t",aux[ k * 16 + l ] ); printf("\n"); - }*/ + } /*bool a = Communicators::MpiCommunicator::IsInitialized(); if( a ) -- GitLab From 7ce9037578f89eb241a22f3df5ce5bbfaa419ca0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matou=C5=A1=20Fencl?= Date: Fri, 30 Nov 2018 15:16:19 +0100 Subject: [PATCH 6/6] MPI ready in tnlDirectEikonal* --- .../hamilton-jacobi/HamiltonJacobiProblem.h | 2 -- .../HamiltonJacobiProblem_impl.h | 20 ++++--------------- .../tnlDirectEikonalMethodsBase_impl.h | 6 +++--- .../tnlFastSweepingMethod2D_impl.h | 20 +++++++++---------- 4 files changed, 17 insertions(+), 31 deletions(-) diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem.h index 7f1bd4193..a41442000 100644 --- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem.h +++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem.h @@ -25,8 +25,6 @@ #include #include -#include - template< typename Mesh, typename DifferentialOperator, typename BoundaryCondition, diff --git a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem_impl.h b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem_impl.h index 9244b1833..3cc638849 100644 --- a/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem_impl.h +++ b/src/TNL/Experimental/Hamilton-Jacobi/Solvers/hamilton-jacobi/HamiltonJacobiProblem_impl.h @@ -123,24 +123,12 @@ setInitialCondition( const Config::ParameterContainer& parameters, { this->bindDofs( mesh, dofs ); const String& initialConditionFile = parameters.getParameter< String >( "initial-condition" ); - if(CommunicatorType::isDistributed()) + if( ! this->solution.boundLoad( initialConditionFile ) ) { - std::cout<<"Nodes Distribution: " << uPointer->getMesh().getDistributedMesh()->printProcessDistr() << std::endl; - if(distributedIOType==Meshes::DistributedMeshes::MpiIO) - Meshes::DistributedMeshes::DistributedGridIO ::load(initialConditionFile, *uPointer ); - if(distributedIOType==Meshes::DistributedMeshes::LocalCopy) - Meshes::DistributedMeshes::DistributedGridIO ::load(initialConditionFile, *uPointer ); - uPointer->template synchronize(); + std::cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." <solution.boundLoad( initialConditionFile ) ) - { - std::cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." < >:: updateBlocks( InterfaceMapType interfaceMap, MeshFunctionType aux, MeshFunctionType helpFunc, - ArrayContainer BlockIterHost, IndexType numThreadsPerBlock/*, Real **sArray*/ ) + ArrayContainer BlockIterHost, int numThreadsPerBlock/*, Real **sArray*/ ) { #pragma omp parallel for schedule( dynamic ) for( IndexType i = 0; i < BlockIterHost.getSize(); i++ ) @@ -267,13 +267,13 @@ updateBlocks( InterfaceMapType interfaceMap, template< typename Real, typename Device, typename Index > -template< IndexType sizeSArray > +template< int sizeSArray > void tnlDirectEikonalMethodsBase< Meshes::Grid< 3, Real, Device, Index > >:: updateBlocks( const InterfaceMapType interfaceMap, const MeshFunctionType aux, MeshFunctionType& helpFunc, - ArrayContainer BlockIterHost, IndexType numThreadsPerBlock/*, Real **sArray*/ ) + ArrayContainer BlockIterHost, int numThreadsPerBlock/*, Real **sArray*/ ) { //#pragma omp parallel for schedule( dynamic ) for( IndexType i = 0; i < BlockIterHost.getSize(); i++ ) 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 27d266341..001f4c69a 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 @@ -85,25 +85,25 @@ solve( const MeshPointer& mesh, #ifdef HAVE_MPI int i = Communicators::MpiCommunicator::GetRank( Communicators::MpiCommunicator::AllGroup ); - printf( "Hello world from rank: %d ", i ); + //printf( "Hello world from rank: %d ", i ); //Communicators::MpiCommunicator::Request r = Communicators::MpiCommunicator::ISend( auxPtr, 0, 0, Communicators::MpiCommunicator::AllGroup ); - /*if( i == 1 ) - for( int k = 0; k < 16*16; k++ ) - aux[ k ] = 10; + if( i == 1 ) + /*for( int k = 0; k < 16*16; k++ ) + aux[ k ] = 10;*/ printf( "1: mesh x: %d\n", mesh->getDimensions().x() ); printf( "1: mesh y: %d\n", mesh->getDimensions().y() ); - aux.save("aux_proc1.tnl");*/ + //aux.save("aux_proc1.tnl"); if( i == 0 ) printf( "0: mesh x: %d\n", mesh->getDimensions().x() ); printf( "0: mesh y: %d\n", mesh->getDimensions().y() ); - aux.save("aux_proc0.tnl"); - for( int k = 0; k < 16*16; k++ ) + //aux.save("aux_proc0.tnl"); + /*for( int k = 0; k < mesh->getDimensions().x()*mesh->getDimensions().y(); k++ ) aux[ k ] = 10; - for( int k = 0; k < 16; k++ ){ - for( int l = 0; l < 16; l++ ) + for( int k = 0; k < mesh->getDimensions().x(); k++ ){ + for( int l = 0; l < mesh->getDimensions().y(); l++ ) printf("%f.2\t",aux[ k * 16 + l ] ); printf("\n"); - } + }*/ /*bool a = Communicators::MpiCommunicator::IsInitialized(); if( a ) -- GitLab