diff --git a/CMakeLists.txt b/CMakeLists.txt index 88f76b88..3ec21350 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -15,7 +15,7 @@ SET( TEST_MAX_PROCS 16 ) # Initialize the project -PROJECT( ${PROJ} ) +PROJECT( ${PROJ} LANGUAGES CXX ) # Prevent users from building in place @@ -25,8 +25,13 @@ ENDIF() # Set the default C++ standard -IF ( NOT CXX_STD ) - SET( CXX_STD 11 ) +SET( CMAKE_CXX_EXTENSIONS OFF ) +IF ( NOT CMAKE_CXX_STANDARD ) + IF ( CXX_STD ) + MESSAGE( WARNING "CXX_STD is obsolete, please set CMAKE_CXX_STANDARD" ) + SET( CMAKE_CXX_STANDARD ${CXX_STD} ) + ENDIF() + SET( CMAKE_CXX_STANDARD 11 ) ENDIF() @@ -113,6 +118,8 @@ ADD_DISTCLEAN( analysis null_timer tests liblbpm-wia.* cpu gpu example common vi # Check for CUDA +CHECK_ENABLE_FLAG( USE_CUDA 0 ) +NULL_USE( CMAKE_CUDA_FLAGS ) IF ( USE_CUDA ) ADD_DEFINITIONS( -DUSE_CUDA ) ENABLE_LANGUAGE( CUDA ) diff --git a/IO/Writer.cpp b/IO/Writer.cpp index 12bf07d0..bb522cf6 100644 --- a/IO/Writer.cpp +++ b/IO/Writer.cpp @@ -121,7 +121,7 @@ static std::vector writeMeshesOrigFormat( const std::vector N = { mesh.nx, mesh.ny, mesh.nz }; - std::array baseindex = { info.ix, info.jy, info.kz }; - const std::string meshname = database.domains[0].name; + auto meshname = database.domains[0].name; silo::writeUniformMesh<3>( fid, meshname, range, N ); silo::write( fid, meshname+"_rankinfo", { mesh.rank, mesh.nprocx, mesh.nprocy, mesh.nprocz } ); for (size_t i=0; i( var.type ); if ( var.precision == IO::DataType::Double ) { silo::writeUniformMeshVariable<3>( fid, meshname, N, var.name, var.data, type ); @@ -336,10 +335,8 @@ static void writeSiloDomainMesh( DBfile *fid, const IO::MeshDataStruct& meshData static IO::MeshDatabase write_domain_silo( DBfile *fid, const std::string& filename, const IO::MeshDataStruct& mesh, int format ) { - const int level = 0; - int rank = MPI_WORLD_RANK(); // Create the MeshDatabase - IO::MeshDatabase database = getDatabase( filename, mesh, format ); + auto database = getDatabase( filename, mesh, format ); if ( database.meshClass=="PointList" ) { writeSiloPointList( fid, mesh, database ); } else if ( database.meshClass=="TriMesh" ) { @@ -429,7 +426,7 @@ static std::vector writeMeshesSilo( sprintf(fullpath,"%s/%s",path.c_str(),filename); auto fid = silo::open( fullpath, silo::CREATE ); for (size_t i=0; i mesh = meshData[i].mesh; + auto mesh = meshData[i].mesh; meshes_written.push_back( write_domain_silo(fid,filename,meshData[i],format) ); } silo::close( fid ); diff --git a/analysis/eikonal.hpp b/analysis/eikonal.cpp similarity index 95% rename from analysis/eikonal.hpp rename to analysis/eikonal.cpp index bd539ea1..82d46246 100644 --- a/analysis/eikonal.hpp +++ b/analysis/eikonal.cpp @@ -1,13 +1,10 @@ -#ifndef Eikonal_HPP_INC -#define Eikonal_HPP_INC - #include "analysis/eikonal.h" #include "analysis/imfilter.h" -inline float minmod(float &a, float &b) +static inline float minmod(float &a, float &b) { float value = a; if ( a*b < 0.0) @@ -18,7 +15,7 @@ inline float minmod(float &a, float &b) } -inline double minmod(double &a, double &b){ +static inline double minmod(double &a, double &b){ double value; @@ -33,9 +30,8 @@ inline double minmod(double &a, double &b){ /****************************************************************** * Solve the eikonal equation * ******************************************************************/ - - -inline double Eikonal(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){ +double Eikonal(DoubleArray &Distance, const char *ID, const Domain &Dm, int timesteps) +{ /* * This routine converts the data in the Distance array to a signed distance @@ -180,7 +176,7 @@ inline double Eikonal(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps return GlobalVar; } -inline float Eikonal3D( Array &Distance, const Array &ID, const Domain &Dm, const int timesteps) +float Eikonal3D( Array &Distance, const Array &ID, const Domain &Dm, const int timesteps) { PROFILE_START("Eikonal3D"); @@ -327,7 +323,7 @@ inline float Eikonal3D( Array &Distance, const Array &ID, const Dom /****************************************************************** * A fast distance calculation * ******************************************************************/ -inline bool CalcDist3DIteration( Array &Distance, const Domain &Dm ) +bool CalcDist3DIteration( Array &Distance, const Domain &Dm ) { const float sq2 = sqrt(2.0f); const float sq3 = sqrt(3.0f); @@ -367,7 +363,7 @@ inline bool CalcDist3DIteration( Array &Distance, const Domain &Dm ) } return changed; } -inline void CalcDist3D( Array &Distance, const Array &ID, const Domain &Dm ) +void CalcDist3D( Array &Distance, const Array &ID, const Domain &Dm ) { PROFILE_START("Calc Distance"); // Initialize the distance to be 0 fore the cells adjacent to the interface @@ -402,7 +398,7 @@ inline void CalcDist3D( Array &Distance, const Array &ID, const Dom /****************************************************************** * A fast distance calculation * ******************************************************************/ -inline void CalcDistMultiLevelHelper( Array &Distance, const Domain &Dm ) +void CalcDistMultiLevelHelper( Array &Distance, const Domain &Dm ) { size_t ratio = 4; std::function&)> coarsen = [ratio]( const Array& data ) @@ -440,7 +436,10 @@ inline void CalcDistMultiLevelHelper( Array &Distance, const Domain &Dm ) // Coarsen Array dist(Nx,Ny,Nz); fillData.copy(Distance,dist); - Domain Dm2(Nx/ratio,Ny/ratio,Nz/ratio,Dm.rank,Dm.nprocx,Dm.nprocy,Dm.nprocz,Dm.Lx,Dm.Ly,Dm.Lz,0); + auto db = Dm.getDatabase()->cloneDatabase(); + auto n = db->getVector( "n" ); + db->putVector( "n", { n[0]/ratio, n[1]/ratio, n[2]/ratio } ); + Domain Dm2(db); Dm2.CommInit(Dm.Comm); fillHalo fillData2(Dm2.Comm,Dm2.rank_info,Nx/ratio,Ny/ratio,Nz/ratio,1,1,1,0,1); auto dist2 = dist.coarsen( {ratio,ratio,ratio}, coarsen ); @@ -483,7 +482,7 @@ inline void CalcDistMultiLevelHelper( Array &Distance, const Domain &Dm ) } } } -inline void CalcDistMultiLevel( Array &Distance, const Array &ID, const Domain &Dm ) +void CalcDistMultiLevel( Array &Distance, const Array &ID, const Domain &Dm ) { PROFILE_START("Calc Distance Multilevel"); int Nx = Dm.Nx-2; @@ -515,5 +514,3 @@ inline void CalcDistMultiLevel( Array &Distance, const Array &ID, c Distance = imfilter::imfilter_separable( Distance, {H,H,H}, BC ); PROFILE_STOP("Calc Distance Multilevel"); } - -#endif diff --git a/analysis/eikonal.h b/analysis/eikonal.h index 6cc0896c..5157e1c5 100644 --- a/analysis/eikonal.h +++ b/analysis/eikonal.h @@ -1,6 +1,7 @@ #ifndef Eikonal_H_INC #define Eikonal_H_INC +#include "common/Domain.h" /*! @@ -21,7 +22,7 @@ * @param[in] timesteps Maximum number of timesteps to process * @return Returns the global variation */ -inline double Eikonal(DoubleArray &Distance, const char *ID, const Domain &Dm, int timesteps); +double Eikonal(DoubleArray &Distance, const char *ID, const Domain &Dm, int timesteps); float Eikonal3D( Array &Distance, const Array &ID, const Domain &Dm, const int timesteps); @@ -52,7 +53,4 @@ void CalcDist3D( Array &Distance, const Array &ID, const Domain &Dm void CalcDistMultiLevel( Array &Distance, const Array &ID, const Domain &Dm ); - -#include "analysis/eikonal.hpp" - #endif diff --git a/cmake/macros.cmake b/cmake/macros.cmake index ddb8db85..374fcd66 100644 --- a/cmake/macros.cmake +++ b/cmake/macros.cmake @@ -35,11 +35,14 @@ IF ( ${CMAKE_BUILD_TYPE} STREQUAL "Release" AND ${CMAKE_VERSION} VERSION_GREATER CHECK_IPO_SUPPORTED(RESULT supported OUTPUT error) IF( supported ) MESSAGE(STATUS "IPO / LTO enabled") + SET( LTO TRUE ) SET( CMAKE_INTERPROCEDURAL_OPTIMIZATION TRUE ) ELSE() - MESSAGE(STATUS "IPO / LTO not supported: <${error}>") + SET( LTO FALSE ) + MESSAGE(STATUS "IPO / LTO not supported") ENDIF() ELSEIF( NOT DEFINED CMAKE_INTERPROCEDURAL_OPTIMIZATION ) + SET( LTO FALSE ) MESSAGE(STATUS "IPO / LTO disabled") ENDIF() @@ -320,7 +323,7 @@ MACRO( INSTALL_PROJ_LIB ) FOREACH ( tmp ${${PROJ}_LIBS} ) SET( tmp_link_list ${tmp_link_list} $ ) ENDFOREACH() - ADD_LIBRARY( ${${PROJ}_LIB} ${tmp_link_list} ) + ADD_LIBRARY( ${${PROJ}_LIB} ${LIB_TYPE} ${tmp_link_list} ) TARGET_LINK_EXTERNAL_LIBRARIES( ${${PROJ}_LIB} LINK_PUBLIC ) INSTALL( TARGETS ${${PROJ}_LIB} DESTINATION ${${PROJ}_INSTALL_DIR}/lib ) ENDMACRO() @@ -408,6 +411,9 @@ MACRO( IDENTIFY_COMPILER ) IF( CMAKE_COMPILER_IS_GNUG77 OR (${CMAKE_Fortran_COMPILER_ID} MATCHES "GNU") ) SET( USING_GFORTRAN TRUE ) MESSAGE("Using gfortran") + IF ( NOT USING_GCC ) + LIST( REMOVE_ITEM CMAKE_Fortran_IMPLICIT_LINK_LIBRARIES gcc ) + ENDIF() ELSEIF ( (${CMAKE_Fortran_COMPILER_ID} MATCHES "Intel") ) SET(USING_IFORT TRUE) MESSAGE("Using ifort") @@ -477,76 +483,6 @@ MACRO( ADD_USER_FLAGS ) ENDMACRO() -# Macro to add user c++ std -MACRO( ADD_CXX_STD ) - IF ( NOT CXX_STD ) - MESSAGE( FATAL_ERROR "The desired c++ standard must be specified: CXX_STD=(98,11,14,NONE)") - ENDIF() - IF ( ${CXX_STD} STREQUAL "NONE" ) - # Do nothing - return() - ELSEIF ( (NOT ${CXX_STD} STREQUAL "98") AND (NOT ${CXX_STD} STREQUAL "11") AND (NOT ${CXX_STD} STREQUAL "14") ) - MESSAGE( FATAL_ERROR "Unknown c++ standard (98,11,14,NONE)" ) - ENDIF() - # Add the flags - SET( CMAKE_CXX_STANDARD ${CXX_STD} ) - MESSAGE( "C++ standard: ${CXX_STD}" ) - SET( CXX_STD_FLAG ) - IF ( USING_GCC ) - # GNU: -std= - IF ( ${CXX_STD} STREQUAL "98" ) - SET( CXX_STD_FLAG -std=c++98 ) - ELSEIF ( ${CXX_STD} STREQUAL "11" ) - SET( CXX_STD_FLAG -std=c++11 ) - ELSEIF ( ${CXX_STD} STREQUAL "14" ) - SET( CXX_STD_FLAG -std=c++1y ) - ELSE() - MESSAGE(FATAL_ERROR "Unknown standard") - ENDIF() - ELSEIF ( USING_MSVC ) - # Microsoft: Does not support this level of control - ELSEIF ( USING_ICC ) - # ICC: -std= - SET( CXX_STD_FLAG -std=c++${CXX_STD} ) - ELSEIF ( USING_CRAY ) - # Cray: Does not seem to support controlling the std? - ELSEIF ( USING_PGCC ) - # PGI: -std= - IF ( ${CXX_STD} STREQUAL "98" ) - SET( CXX_STD_FLAG --c++0x ) - ELSEIF ( ${CXX_STD} STREQUAL "11" ) - SET( CXX_STD_FLAG --c++11 ) - ELSEIF ( ${CXX_STD} STREQUAL "14" ) - MESSAGE( FATAL_ERROR "C++14 features are not availible yet for PGI" ) - ELSE() - MESSAGE(FATAL_ERROR "Unknown standard") - ENDIF() - ELSEIF ( USING_CLANG ) - # Clang: -std= - IF ( ( ${CXX_STD} STREQUAL "98") OR ( ${CXX_STD} STREQUAL "11" ) ) - SET( CXX_STD_FLAG -std=c++${CXX_STD} ) - ELSEIF ( ${CXX_STD} STREQUAL "14" ) - SET( CXX_STD_FLAG -std=c++1y ) - ELSE() - MESSAGE(FATAL_ERROR "Unknown standard") - ENDIF() - ELSEIF ( USING_XL ) - # XL: -std= - IF ( ( ${CXX_STD} STREQUAL "98") OR ( ${CXX_STD} STREQUAL "11" ) ) - SET( CXX_STD_FLAG -std=c++${CXX_STD} ) - ELSEIF ( ${CXX_STD} STREQUAL "14" ) - SET( CXX_STD_FLAG -std=c++1y ) - ELSE() - MESSAGE(FATAL_ERROR "Unknown standard") - ENDIF() - ELSEIF ( USING_DEFAULT ) - # Default: do nothing - ENDIF() - ADD_DEFINITIONS( -DCXX_STD=${CXX_STD} ) - SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${CXX_STD_FLAG}") -ENDMACRO() - - # Macro to set the compile/link flags MACRO( SET_COMPILER_FLAGS ) # Initilaize the compiler @@ -585,8 +521,6 @@ MACRO( SET_COMPILER_FLAGS ) SET( CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}" CACHE STRING "Release flags" FORCE) # Add the user flags ADD_USER_FLAGS() - # Add the c++ standard flags - ADD_CXX_STD() # Set the warnings to use SET_WARNINGS() # Test the compile flags @@ -1098,71 +1032,19 @@ ENDMACRO() # Add a matlab mex file -FUNCTION( ADD_MATLAB_MEX MEXFILE ) - # Set the MEX compiler and default link flags - IF ( USING_MSVC ) - SET( MEX_LDFLAGS ${MEX_LDFLAGS} -L${CMAKE_CURRENT_BINARY_DIR}/.. - -L${CMAKE_CURRENT_BINARY_DIR}/../Debug -L${CMAKE_CURRENT_BINARY_DIR} ) - ENDIF() - SET( MEX_LDFLAGS ${MEX_LDFLAGS} -L${CMAKE_CURRENT_BINARY_DIR}/.. -L${CMAKE_CURRENT_BINARY_DIR}/../.. ${SYSTEM_LDFLAGS} ) - SET( MEX_LIBS ${MEX_LIBS} ${SYSTEM_LIBS} ) - IF ( NOT USING_MSVC ) - FOREACH( rpath ${CMAKE_INSTALL_RPATH} ) - SET( MEX_LDFLAGS ${MEX_LDFLAGS} "-Wl,-rpath,${rpath},--no-undefined" ) - ENDFOREACH() - ENDIF() - IF ( "${CMAKE_BUILD_TYPE}" STREQUAL "Debug" ) - SET( MEX_FLAGS ${MEX_FLAGS} -g ) - ELSEIF ( "${CMAKE_BUILD_TYPE}" STREQUAL "Release" ) - SET( MEX_FLAGS ${MEX_FLAGS} -O ) - ELSE() - MESSAGE( FATAL_ERROR "Unknown build type: ${CMAKE_BUILD_TYPE}" ) - ENDIF() - # Create the mex comamnd - STRING(REGEX REPLACE "[.]cpp" "" TARGET ${MEXFILE}) - STRING(REGEX REPLACE "[.]c" "" TARGET ${TARGET}) - FILE(MAKE_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}/temp/${TARGET}" ) - IF ( USING_MSVC ) - SET( MEX_BAT_FILE "${CMAKE_CURRENT_BINARY_DIR}/temp/${TARGET}/compile_mex.bat" ) - FILE( WRITE "${MEX_BAT_FILE}" "${MEX_EXE} \"${CMAKE_CURRENT_SOURCE_DIR}/${MEXFILE}\" " ) - APPEND_LIST( "${MEX_BAT_FILE}" "${MEX_FLAGS}" "\"" "\" " ) - APPEND_LIST( "${MEX_BAT_FILE}" "${COMPFLAGS}" "\"" "\" " ) - APPEND_LIST( "${MEX_BAT_FILE}" "${MEX_INCLUDE}" "\"" "\" " ) - APPEND_LIST( "${MEX_BAT_FILE}" "${MEX_LDFLAGS}" "\"" "\" " ) - APPEND_LIST( "${MEX_BAT_FILE}" "${MEX_LIBS}" "\"" "\" " ) - APPEND_LIST( "${MEX_BAT_FILE}" "${COVERAGE_MATLAB_LIBS}" "\"" "\" " ) - SET( MEX_COMMAND "${MEX_BAT_FILE}" ) - ELSE() - STRING(REPLACE " " ";" MEX_CFLAGS "$$CFLAGS ${CMAKE_C_FLAGS}") - STRING(REPLACE " " ";" MEX_CXXFLAGS "$$CXXFLAGS ${CMAKE_CXX_FLAGS}") - SET( MEX_COMMAND ${MEX_EXE} ${CMAKE_CURRENT_SOURCE_DIR}/${MEXFILE} - ${MEX_FLAGS} ${MEX_INCLUDE} - CFLAGS="${MEX_CFLAGS}" - CXXFLAGS="${MEX_CXXFLAGS}" - LDFLAGS="${MEX_LDFLAGS}" - ${MEX_LIBS} ${COVERAGE_MATLAB_LIBS} - ) - ENDIF() - ADD_CUSTOM_COMMAND( - OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/${TARGET}.${MEX_EXTENSION} - COMMAND ${MEX_COMMAND} - COMMAND ${CMAKE_COMMAND} -E copy "${CMAKE_CURRENT_BINARY_DIR}/temp/${TARGET}/${TARGET}.${MEX_EXTENSION}" "${CMAKE_CURRENT_BINARY_DIR}" - COMMAND ${CMAKE_COMMAND} -E copy "${CMAKE_CURRENT_BINARY_DIR}/temp/${TARGET}/${TARGET}.${MEX_EXTENSION}" "${${PROJ}_INSTALL_DIR}/mex" - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/temp/${TARGET} - DEPENDS ${${PROJ}_LIBS} ${MATLAB_TARGET} ${CMAKE_CURRENT_SOURCE_DIR}/${MEXFILE} +FUNCTION( ADD_MATLAB_MEX SOURCE ) + STRING( REGEX REPLACE "[.]cpp" "" TARGET ${SOURCE} ) + STRING( REGEX REPLACE "[.]c" "" TARGET ${TARGET} ) + MATLAB_ADD_MEX( + NAME ${TARGET} + SRC ${SOURCE} ) - ADD_CUSTOM_TARGET( ${TARGET} - DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/${TARGET}.${MEX_EXTENSION} - DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/${MEXFILE} - ) - IF ( MATLAB_TARGET ) - ADD_DEPENDENCIES( ${TARGET} ${MATLAB_TARGET} ) - ENDIF() - IF ( ${PROJ}_LIBS ) - ADD_DEPENDENCIES( ${TARGET} ${${PROJ}_LIBS} ) - ENDIF() + TARGET_LINK_LIBRARIES( ${TARGET} ${MATLAB_TARGET} ) + ADD_PROJ_EXE_DEP( ${TARGET} ) ADD_DEPENDENCIES( mex ${TARGET} ) - SET( MEX_FILES2 ${MEX_FILES} "${${PROJ}_INSTALL_DIR}/mex/${TARGET}.${MEX_EXTENSION}" ) + INSTALL( TARGETS ${TARGET} DESTINATION ${${PROJ}_INSTALL_DIR}/mex ) + ADD_DEPENDENCIES( mex ${TARGET} ) + SET( MEX_FILES2 ${MEX_FILES} "${${PROJ}_INSTALL_DIR}/mex/${TARGET}.${Matlab_MEX_EXTENSION}" ) LIST( REMOVE_DUPLICATES MEX_FILES2 ) SET( MEX_FILES ${MEX_FILES2} CACHE INTERNAL "" ) ENDFUNCTION() @@ -1178,7 +1060,7 @@ MACRO( ADD_MATLAB_TEST EXEFILE ${ARGN} ) IF ( USING_MSVC ) SET( MATLAB_OPTIONS "-logfile" "log_${EXEFILE}" ) ENDIF() - SET( MATLAB_COMMAND "try, ${EXEFILE}, catch ME, ME, clear all global, exit(1), end, disp('ALL TESTS PASSED'); exit(0)" ) + SET( MATLAB_COMMAND "try, ${EXEFILE}, catch ME, disp(getReport(ME)), clear all global, exit(1), end, disp('ALL TESTS PASSED'); exit(0)" ) SET( MATLAB_DEBUGGER_OPTIONS ) IF ( MATLAB_DEBUGGER ) SET( MATLAB_DEBUGGER_OPTIONS -D${MATLAB_DEBUGGER} ) @@ -1195,7 +1077,6 @@ FUNCTION( CREATE_MATLAB_WRAPPER ) SET( tmp_libs ${MEX_LIBCXX} ${MEX_FILES} ) STRING(REGEX REPLACE ";" ":" tmp_libs "${tmp_libs}") STRING(REGEX REPLACE ";" ":" tmp_path "${MATLABPATH}") - FIND_PROGRAM( MATLAB_EXE2 NAME matlab PATHS "${MATLAB_DIRECTORY}/bin" NO_DEFAULT_PATH ) IF ( USING_MSVC ) # Create a matlab wrapper for windows SET( MATLAB_GUI "${CMAKE_CURRENT_BINARY_DIR}/tmp/matlab-gui.bat" ) @@ -1212,8 +1093,8 @@ FUNCTION( CREATE_MATLAB_WRAPPER ) SET( MATLAB_GUI "${CMAKE_CURRENT_BINARY_DIR}/tmp/matlab-gui" ) SET( MATLAB_CMD "${CMAKE_CURRENT_BINARY_DIR}/tmp/matlab-cmd" ) SET( MATLAB_INSTALL_CMD "matlab-cmd" ) - FILE( WRITE "${MATLAB_GUI}" "LD_PRELOAD=\"${tmp_libs}\" MKL_NUM_THREADS=1 MATLABPATH=\"${tmp_path}\" \"${MATLAB_EXE2}\" -singleCompThread -nosplash \"$@\"\n") - FILE( WRITE "${MATLAB_CMD}" "LD_PRELOAD=\"${tmp_libs}\" MKL_NUM_THREADS=1 MATLABPATH=\"${tmp_path}\" \"${MATLAB_EXE2}\" -singleCompThread -nosplash -nodisplay -nojvm \"$@\"\n") + FILE( WRITE "${MATLAB_GUI}" "LD_PRELOAD=\"${tmp_libs}\" MKL_NUM_THREADS=1 MATLABPATH=\"${tmp_path}\" \"${Matlab_MAIN_PROGRAM}\" -singleCompThread -nosplash \"$@\"\n") + FILE( WRITE "${MATLAB_CMD}" "LD_PRELOAD=\"${tmp_libs}\" MKL_NUM_THREADS=1 MATLABPATH=\"${tmp_path}\" \"${Matlab_MAIN_PROGRAM}\" -singleCompThread -nosplash -nodisplay -nojvm \"$@\"\n") ENDIF() FILE( COPY "${MATLAB_GUI}" DESTINATION "${${PROJ}_INSTALL_DIR}/mex" FILE_PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE ) @@ -1294,6 +1175,7 @@ MACRO( ADD_DISTCLEAN ${ARGN} ) install_manifest.txt test matlab + Matlab mex tmp #tmp# diff --git a/common/Database.cpp b/common/Database.cpp new file mode 100644 index 00000000..d318cc03 --- /dev/null +++ b/common/Database.cpp @@ -0,0 +1,597 @@ +#include "common/Database.h" +#include "common/Utilities.h" + +#include +#include +#include +#include +#include +#include + + +/******************************************************************** + * Constructors/destructor * + ********************************************************************/ +Database::Database() = default; +Database::~Database() = default; +Database::Database( const Database& rhs ) : KeyData( rhs ) +{ + d_data.clear(); + for ( const auto& tmp : rhs.d_data ) + putData( tmp.first, tmp.second->clone() ); +} +Database& Database::operator=( const Database& rhs ) +{ + if ( this == &rhs ) + return *this; + d_data.clear(); + for ( const auto& tmp : rhs.d_data ) + putData( tmp.first, tmp.second->clone() ); + return *this; +} +Database::Database( Database&& rhs ) { std::swap( d_data, rhs.d_data ); } +Database& Database::operator=( Database&& rhs ) +{ + if ( this != &rhs ) + std::swap( d_data, rhs.d_data ); + return *this; +} + + +/******************************************************************** + * Clone the database * + ********************************************************************/ +std::shared_ptr Database::clone() const { return cloneDatabase(); } +std::shared_ptr Database::cloneDatabase() const +{ + auto db = std::make_shared(); + for ( const auto& tmp : d_data ) + db->putData( tmp.first, tmp.second->clone() ); + return db; +} + + +/******************************************************************** + * Get the data object * + ********************************************************************/ +bool Database::keyExists( const std::string& key ) const +{ + return d_data.find( key ) != d_data.end(); +} +std::shared_ptr Database::getData( const std::string& key ) +{ + auto it = d_data.find( key ); + if ( it == d_data.end() ) { + char msg[1000]; + sprintf( msg, "Variable %s was not found in database", key.c_str() ); + ERROR( msg ); + } + return it->second; +} +std::shared_ptr Database::getData( const std::string& key ) const +{ + return const_cast( this )->getData( key ); +} +bool Database::isDatabase( const std::string& key ) const +{ + auto ptr = getData( key ); + auto ptr2 = std::dynamic_pointer_cast( ptr ); + return ptr2 != nullptr; +} +std::shared_ptr Database::getDatabase( const std::string& key ) +{ + std::shared_ptr ptr = getData( key ); + std::shared_ptr ptr2 = std::dynamic_pointer_cast( ptr ); + if ( ptr2 == nullptr ) { + char msg[1000]; + sprintf( msg, "Variable %s is not a database", key.c_str() ); + ERROR( msg ); + } + return ptr2; +} +std::shared_ptr Database::getDatabase( const std::string& key ) const +{ + return const_cast( this )->getDatabase( key ); +} +std::vector Database::getAllKeys() const +{ + std::vector keys; + keys.reserve( d_data.size() ); + for ( const auto& it : d_data ) + keys.push_back( it.first ); + return keys; +} +void Database::putDatabase( const std::string& key, std::shared_ptr db ) +{ + d_data[key] = std::move( db ); +} +void Database::putData( const std::string& key, std::shared_ptr data ) +{ + d_data[key] = std::move( data ); +} + + +/******************************************************************** + * Is the data of the given type * + ********************************************************************/ +template<> +bool Database::isType( const std::string& key ) const +{ + auto type = getData( key )->type(); + return type == "double"; +} +template<> +bool Database::isType( const std::string& key ) const +{ + auto type = getData( key )->type(); + return type == "double"; +} +template<> +bool Database::isType( const std::string& key ) const +{ + bool pass = true; + auto type = getData( key )->type(); + if ( type == "double" ) { + auto data = getVector( key ); + for ( auto tmp : data ) + pass = pass && static_cast( static_cast( tmp ) ) == tmp; + } else { + pass = false; + } + return pass; +} +template<> +bool Database::isType( const std::string& key ) const +{ + auto type = getData( key )->type(); + return type == "string"; +} +template<> +bool Database::isType( const std::string& key ) const +{ + auto type = getData( key )->type(); + return type == "bool"; +} + + +/******************************************************************** + * Get a vector * + ********************************************************************/ +template<> +std::vector Database::getVector( + const std::string& key, const Units& ) const +{ + std::shared_ptr ptr = getData( key ); + if ( std::dynamic_pointer_cast( ptr ) ) + return std::vector(); + const auto* ptr2 = dynamic_cast( ptr.get() ); + if ( ptr2 == nullptr ) { + ERROR( "Key '" + key + "' is not a string" ); + } + return ptr2->d_data; +} +template<> +std::vector Database::getVector( const std::string& key, const Units& ) const +{ + std::shared_ptr ptr = getData( key ); + if ( std::dynamic_pointer_cast( ptr ) ) + return std::vector(); + const auto* ptr2 = dynamic_cast( ptr.get() ); + if ( ptr2 == nullptr ) { + ERROR( "Key '" + key + "' is not a bool" ); + } + return ptr2->d_data; +} +template +std::vector Database::getVector( const std::string& key, const Units& unit ) const +{ + std::shared_ptr ptr = getData( key ); + if ( std::dynamic_pointer_cast( ptr ) ) + return std::vector(); + std::vector data; + if ( std::dynamic_pointer_cast( ptr ) ) { + const auto* ptr2 = dynamic_cast( ptr.get() ); + const std::vector& data2 = ptr2->d_data; + double factor = 1; + if ( !unit.isNull() ) { + INSIST( !ptr2->d_unit.isNull(), "Field " + key + " must have units" ); + factor = ptr2->d_unit.convert( unit ); + INSIST( factor != 0, "Unit conversion failed" ); + } + data.resize( data2.size() ); + for ( size_t i = 0; i < data2.size(); i++ ) + data[i] = static_cast( factor * data2[i] ); + } else if ( std::dynamic_pointer_cast( ptr ) ) { + ERROR( "Converting std::string to another type" ); + } else if ( std::dynamic_pointer_cast( ptr ) ) { + ERROR( "Converting std::bool to another type" ); + } else { + ERROR( "Unable to convert data format" ); + } + return data; +} + + +/******************************************************************** + * Put a vector * + ********************************************************************/ +template<> +void Database::putVector( + const std::string& key, const std::vector& data, const Units& ) +{ + std::shared_ptr ptr( new KeyDataString() ); + ptr->d_data = data; + d_data[key] = ptr; +} +template<> +void Database::putVector( + const std::string& key, const std::vector& data, const Units& ) +{ + std::shared_ptr ptr( new KeyDataBool() ); + ptr->d_data = data; + d_data[key] = ptr; +} +template +void Database::putVector( const std::string& key, const std::vector& data, const Units& unit ) +{ + std::shared_ptr ptr( new KeyDataDouble() ); + ptr->d_unit = unit; + ptr->d_data.resize( data.size() ); + for ( size_t i = 0; i < data.size(); i++ ) + ptr->d_data[i] = static_cast( data[i] ); + d_data[key] = ptr; +} + + +/******************************************************************** + * Print the database * + ********************************************************************/ +void Database::print( std::ostream& os, const std::string& indent ) const +{ + for ( const auto& it : d_data ) { + os << indent << it.first; + if ( dynamic_cast( it.second.get() ) ) { + const auto* db = dynamic_cast( it.second.get() ); + os << " {\n"; + db->print( os, indent + " " ); + os << indent << "}\n"; + } else { + os << " = "; + it.second->print( os, "" ); + } + } +} +std::string Database::print( const std::string& indent ) const +{ + std::stringstream ss; + print( ss, indent ); + return ss.str(); +} + + +/******************************************************************** + * Read input database file * + ********************************************************************/ +Database::Database( const std::string& filename ) +{ + // Read the input file into memory + FILE* fid = fopen( filename.c_str(), "rb" ); + if ( fid == nullptr ) + ERROR( "Error opening file " + filename ); + fseek( fid, 0, SEEK_END ); + size_t bytes = ftell( fid ); + rewind( fid ); + auto* buffer = new char[bytes + 4]; + size_t result = fread( buffer, 1, bytes, fid ); + fclose( fid ); + if ( result != bytes ) + ERROR( "Error reading file " + filename ); + buffer[bytes + 0] = '\n'; + buffer[bytes + 1] = '}'; + buffer[bytes + 2] = '\n'; + buffer[bytes + 3] = 0; + // Create the database entries + loadDatabase( buffer, *this ); + // Free temporary memory + delete[] buffer; +} +std::shared_ptr Database::createFromString( const std::string& data ) +{ + std::shared_ptr db( new Database() ); + auto* buffer = new char[data.size() + 4]; + memcpy( buffer, data.data(), data.size() ); + buffer[data.size() + 0] = '\n'; + buffer[data.size() + 1] = '}'; + buffer[data.size() + 2] = '\n'; + buffer[data.size() + 3] = 0; + loadDatabase( buffer, *db ); + delete[] buffer; + return db; +} +enum class token_type { + newline, + line_comment, + block_start, + block_stop, + quote, + equal, + bracket, + end_bracket, + end +}; +inline size_t length( token_type type ) +{ + size_t len = 0; + if ( type == token_type::newline || type == token_type::quote || type == token_type::equal || + type == token_type::bracket || type == token_type::end_bracket || + type == token_type::end ) { + len = 1; + } else if ( type == token_type::line_comment || type == token_type::block_start || + type == token_type::block_stop ) { + len = 2; + } + return len; +} +inline std::tuple find_next_token( const char* buffer ) +{ + size_t i = 0; + while ( true ) { + if ( buffer[i] == '\n' || buffer[i] == '\r' ) { + return std::pair( i + 1, token_type::newline ); + } else if ( buffer[i] == 0 ) { + return std::pair( i + 1, token_type::end ); + } else if ( buffer[i] == '"' ) { + return std::pair( i + 1, token_type::quote ); + } else if ( buffer[i] == '=' ) { + return std::pair( i + 1, token_type::equal ); + } else if ( buffer[i] == '{' ) { + return std::pair( i + 1, token_type::bracket ); + } else if ( buffer[i] == '}' ) { + return std::pair( i + 1, token_type::end_bracket ); + } else if ( buffer[i] == '/' ) { + if ( buffer[i + 1] == '/' ) { + return std::pair( i + 2, token_type::line_comment ); + } else if ( buffer[i + 1] == '*' ) { + return std::pair( i + 2, token_type::block_start ); + } + } else if ( buffer[i] == '*' ) { + if ( buffer[i + 1] == '/' ) + return std::pair( i + 2, token_type::block_stop ); + } + i++; + } + return std::pair( 0, token_type::end ); +} +inline std::string deblank( const std::string& str ) +{ + size_t i1 = 0xFFFFFFF, i2 = 0; + for ( size_t i = 0; i < str.size(); i++ ) { + if ( str[i] != ' ' ) { + i1 = std::min( i1, i ); + i2 = std::max( i2, i ); + } + } + return i1 <= i2 ? str.substr( i1, i2 - i1 + 1 ) : std::string(); +} +size_t skip_comment( const char* buffer ) +{ + auto tmp = find_next_token( buffer ); + const token_type end_comment = ( std::get<1>( tmp ) == token_type::line_comment ) ? + token_type::newline : + token_type::block_stop; + size_t pos = 0; + while ( std::get<1>( tmp ) != end_comment ) { + if ( std::get<1>( tmp ) == token_type::end ) + ERROR( "Encountered end of file before block comment end" ); + pos += std::get<0>( tmp ); + tmp = find_next_token( &buffer[pos] ); + } + pos += std::get<0>( tmp ); + return pos; +} +inline std::string lower( const std::string& str ) +{ + std::string tmp( str ); + std::transform( tmp.begin(), tmp.end(), tmp.begin(), ::tolower ); + return tmp; +} +static std::tuple> read_value( + const char* buffer, const std::string& key ) +{ + // Get the value as a std::string + size_t pos = 0; + token_type type = token_type::end; + std::tie( pos, type ) = find_next_token( &buffer[pos] ); + size_t len = pos - length( type ); + while ( type != token_type::newline ) { + if ( type == token_type::quote ) { + size_t i = 0; + std::tie( i, type ) = find_next_token( &buffer[pos] ); + pos += i; + while ( type != token_type::quote ) { + ASSERT( type != token_type::end ); + std::tie( i, type ) = find_next_token( &buffer[pos] ); + pos += i; + } + } else if ( type == token_type::line_comment || type == token_type::block_start ) { + len = pos - length( type ); + pos += skip_comment( &buffer[pos - length( type )] ) - length( type ); + break; + } + size_t i = 0; + std::tie( i, type ) = find_next_token( &buffer[pos] ); + pos += i; + len = pos - length( type ); + } + const std::string value = deblank( std::string( buffer, len ) ); + // Split the value to an array of values + std::vector values; + size_t i0 = 0, i = 0, count = 0; + for ( ; i < value.size(); i++ ) { + if ( value[i] == '"' ) { + count++; + } else if ( value[i] == ',' && count % 2 == 0 ) { + values.push_back( deblank( value.substr( i0, i - i0 ) ) ); + i0 = i + 1; + } + } + values.push_back( deblank( value.substr( i0 ) ) ); + // Convert the string value to the database value + std::shared_ptr data; + if ( value.empty() ) { + data.reset( new EmptyKeyData() ); + } else if ( value.find( '"' ) != std::string::npos ) { + auto* data2 = new KeyDataString(); + data.reset( data2 ); + data2->d_data.resize( values.size() ); + for ( size_t i = 0; i < values.size(); i++ ) { + ASSERT( values[i].size() >= 2 ); + ASSERT( values[i][0] == '"' && values[i][values[i].size() - 1] == '"' ); + data2->d_data[i] = values[i].substr( 1, values[i].size() - 2 ); + } + } else if ( lower( value ) == "true" || lower( value ) == "false" ) { + auto* data2 = new KeyDataBool(); + data.reset( data2 ); + data2->d_data.resize( values.size() ); + for ( size_t i = 0; i < values.size(); i++ ) { + ASSERT( values[i].size() >= 2 ); + if ( lower( values[i] ) != "true" && lower( values[i] ) != "false" ) + ERROR( "Error converting " + key + " to logical array" ); + data2->d_data[i] = lower( values[i] ) == "true"; + } + } else { // if ( value.find('.')!=std::string::npos || value.find('e')!=std::string::npos ) { + auto* data2 = new KeyDataDouble(); + data.reset( data2 ); + data2->d_data.resize( values.size(), 0 ); + for ( size_t i = 0; i < values.size(); i++ ) { + Units unit; + std::tie( data2->d_data[i], unit ) = KeyDataDouble::read( values[i] ); + if ( !unit.isNull() ) + data2->d_unit = unit; + } + //} else { + // ERROR("Unable to determine data type: "+value); + } + return std::tuple>( pos, data ); +} +size_t Database::loadDatabase( const char* buffer, Database& db ) +{ + size_t pos = 0; + while ( true ) { + size_t i; + token_type type; + std::tie( i, type ) = find_next_token( &buffer[pos] ); + const std::string key = + deblank( std::string( &buffer[pos], std::max( i - length( type ), 1 ) - 1 ) ); + if ( type == token_type::line_comment || type == token_type::block_start ) { + // Comment + INSIST( key.empty(), "Key should be empty: " + key ); + pos += skip_comment( &buffer[pos] ); + } else if ( type == token_type::newline ) { + INSIST( key.empty(), "Key should be empty: " + key ); + pos += i; + } else if ( type == token_type::equal ) { + // Reading key/value pair + ASSERT( !key.empty() ); + pos += i; + std::shared_ptr data; + std::tie( i, data ) = read_value( &buffer[pos], key ); + ASSERT( data.get() != nullptr ); + db.d_data[key] = data; + pos += i; + } else if ( type == token_type::bracket ) { + // Read database + ASSERT( !key.empty() ); + pos += i; + std::shared_ptr database( new Database() ); + pos += loadDatabase( &buffer[pos], *database ); + db.d_data[key] = database; + } else if ( type == token_type::end_bracket ) { + // Finished with the database + pos += i; + break; + } else { + ERROR( "Error loading data" ); + } + } + return pos; +} + + +/******************************************************************** + * Data type helper functions * + ********************************************************************/ +void KeyDataDouble::print( std::ostream& os, const std::string& indent ) const +{ + os << indent; + for ( size_t i = 0; i < d_data.size(); i++ ) { + if ( i > 0 ) + os << ", "; + if ( d_data[i] != d_data[i] ) { + os << "nan"; + } else if ( d_data[i] == std::numeric_limits::infinity() ) { + os << "inf"; + } else if ( d_data[i] == -std::numeric_limits::infinity() ) { + os << "-inf"; + } else { + os << std::setprecision( 12 ) << d_data[i]; + } + } + if ( !d_unit.isNull() ) + os << " " << d_unit.str(); + os << std::endl; +} +std::tuple KeyDataDouble::read( const std::string& str ) +{ + std::string tmp = deblank( str ); + size_t index = tmp.find( " " ); + if ( index != std::string::npos ) { + return std::make_tuple( + readValue( tmp.substr( 0, index ) ), Units( tmp.substr( index + 1 ) ) ); + } else { + return std::make_tuple( readValue( tmp ), Units() ); + } +} +double KeyDataDouble::readValue( const std::string& str ) +{ + const std::string tmp = lower( str ); + double data = 0; + if ( tmp == "inf" || tmp == "infinity" ) { + data = std::numeric_limits::infinity(); + } else if ( tmp == "-inf" || tmp == "-infinity" ) { + data = -std::numeric_limits::infinity(); + } else if ( tmp == "nan" ) { + data = std::numeric_limits::quiet_NaN(); + } else if ( tmp.find( '/' ) != std::string::npos ) { + ERROR( "Error reading value" ); + } else { + char* pos = nullptr; + data = strtod( tmp.c_str(), &pos ); + if ( static_cast( pos - tmp.c_str() ) == tmp.size() + 1 ) + ERROR( "Error reading value" ); + } + return data; +} + + +/******************************************************************** + * Instantiations * + ********************************************************************/ +template std::vector Database::getVector( const std::string&, const Units& ) const; +template std::vector Database::getVector( const std::string&, const Units& ) const; +template std::vector Database::getVector( const std::string&, const Units& ) const; +template std::vector Database::getVector( const std::string&, const Units& ) const; +template std::vector Database::getVector( const std::string&, const Units& ) const; +template void Database::putVector( + const std::string&, const std::vector&, const Units& ); +template void Database::putVector( const std::string&, const std::vector&, const Units& ); +template void Database::putVector( + const std::string&, const std::vector&, const Units& ); +template void Database::putVector( + const std::string&, const std::vector&, const Units& ); +template void Database::putVector( + const std::string&, const std::vector&, const Units& ); +template bool Database::isType( const std::string& ) const; +template bool Database::isType( const std::string& ) const; +template bool Database::isType( const std::string& ) const; +template bool Database::isType( const std::string& ) const; diff --git a/common/Database.h b/common/Database.h new file mode 100644 index 00000000..4bf0bf61 --- /dev/null +++ b/common/Database.h @@ -0,0 +1,288 @@ +#ifndef included_Database +#define included_Database + +#include +#include +#include +#include +#include + +#include "common/Units.h" + + +//! Base class to hold data of a given type +class KeyData +{ +protected: + //! Empty constructor + KeyData() {} + +public: + //! Destructor + virtual ~KeyData() {} + //! Copy the data + virtual std::shared_ptr clone() const = 0; + //! Print the data to a stream + virtual void print( std::ostream& os, const std::string& indent = "" ) const = 0; + //! Return the native data type + virtual std::string type() const = 0; + +protected: + KeyData( const KeyData& ) {} + KeyData& operator=( const KeyData& ); +}; + + +//! Class to a database +class Database : public KeyData +{ +public: + //! Empty constructor + Database(); + + /** + * Open an database file. + * @param filename Name of input file to open + */ + explicit Database( const std::string& filename ); + + /** + * Create database from string + * @param data String containing the database data + */ + static std::shared_ptr createFromString( const std::string& data ); + + //! Copy constructor + Database( const Database& ); + + //! Assignment operator + Database& operator=( const Database& ); + + //! Move constructor + Database( Database&& rhs ); + + //! Move assignment operator + Database& operator=( Database&& rhs ); + + //! Destructor + virtual ~Database(); + + //! Copy the data + virtual std::shared_ptr clone() const override; + + //! Copy the data + std::shared_ptr cloneDatabase() const; + + + /** + * Return true if the specified key exists in the database and false + * otherwise. + * @param[in] key Key name to lookup. + */ + bool keyExists( const std::string& key ) const; + + + /** + * Return all keys in the database. + */ + std::vector getAllKeys() const; + + + //! Return the number of entries in the database + size_t size() const { return d_data.size(); } + + + /** + * Get the scalar entry from the database with the specified key + * name. If the specified key does not exist in the database or + * is not a scalar of the given type, then an error message is printed and + * the program exits. + * + * @param[in] key Key name in database. + * @param[in] unit Desired units + */ + template + inline TYPE getScalar( const std::string& key, const Units& unit = Units() ) const; + + + /// @copydoc Database::getScalar(const std::string&,const Units&) const + template + inline TYPE getScalar( const std::string& key, const std::string& unit ) const; + + + /** + * Get the scalar entry from the database with the specified key + * name. If the specified key does not exist in the database the + * the default value will be printed + * + * @param[in] key Key name in database + * @param[in] value Default value + * @param[in] unit Desired units + */ + template + inline TYPE getWithDefault( + const std::string& key, const TYPE& value, const Units& unit = Units() ) const; + + + /// @copydoc Database::getWithDefault(const std::string&,const TYPE&,const Units&) const + template + inline TYPE getWithDefault( + const std::string& key, const TYPE& value, const std::string& unit ) const; + + + /** + * Put the scalar entry into the database with the specified key name. + * @param key Key name in database. + * @param value Value to store + * @param unit Desired units + */ + template + inline void putScalar( const std::string& key, const TYPE& value, const Units& unit = Units() ); + + + /** + * Put the scalar entry into the database with the specified key name. + * @param key Key name in database. + * @param value Value to store + * @param unit Desired units + */ + template + inline void putScalar( const std::string& key, const TYPE& value, const std::string& unit ); + + + /** + * Get the vector entries from the database with the specified key + * name. If the specified key does not exist in the database or + * is not of the given type, then an error message is printed and + * the program exits. + * + * @param key Key name in database. + * @param unit Desired units + */ + template + std::vector getVector( const std::string& key, const Units& unit = Units() ) const; + + + /// @copydoc Database::getVector(const std::string&,const Units&) const + template + inline std::vector getVector( const std::string& key, const std::string& unit ) const; + + + /** + * Put the vector entries into the database with the specified key + * name. If the specified key does not exist in the database or + * is not of the given type, then an error message is printed and + * the program exits. + * + * @param key Key name in database. + * @param data Data to store + * @param unit Desired units + */ + template + void putVector( + const std::string& key, const std::vector& data, const Units& unit = Units() ); + + + /// @copydoc Database::putVector(const std::string&,const std::vector&,const Units&) + template + inline void putVector( + const std::string& key, const std::vector& data, const std::string& unit ); + + + /** + * Get the data for a key in the database. If the specified key + * does not exist in the database an error message is printed and + * the program exits. + * + * @param key Key name in database. + */ + std::shared_ptr getData( const std::string& key ); + + /** + * Get the data for a key in the database. If the specified key + * does not exist in the database an error message is printed and + * the program exits. + * + * @param key Key name in database. + */ + std::shared_ptr getData( const std::string& key ) const; + + + /** + * Put the data for a key in the database. + * + * @param key Key name in database. + * @param data Data to store + */ + void putData( const std::string& key, std::shared_ptr data ); + + + // Check if the key is a database object + bool isDatabase( const std::string& key ) const; + + + // Check if the entry can be stored as the given type + template + bool isType( const std::string& key ) const; + + + /** + * Get the database for a key in the database. If the specified key + * does not exist in the database an error message is printed and + * the program exits. + * + * @param key Key name in database. + */ + std::shared_ptr getDatabase( const std::string& key ); + + /** + * Get the database for a key in the database. If the specified key + * does not exist in the database an error message is printed and + * the program exits. + * + * @param key Key name in database. + */ + std::shared_ptr getDatabase( const std::string& key ) const; + + + /** + * Get the database for a key in the database. If the specified key + * does not exist in the database an error message is printed and + * the program exits. + * + * @param key Key name in database. + * @param db Database to store + */ + void putDatabase( const std::string& key, std::shared_ptr db ); + + + /** + * Print the data to a stream + * @param os Output stream + * @param indent Indenting to use before each line + */ + virtual void print( std::ostream& os, const std::string& indent = "" ) const override; + + + //! Print the type + virtual std::string type() const override { return "database"; } + + + /** + * Print the data to a string + * @return Output string + */ + std::string print( const std::string& indent = "" ) const; + + +protected: + std::map> d_data; + + // Function to load a database from a buffer + static size_t loadDatabase( const char* buffer, Database& db ); +}; + + +#include "common/Database.hpp" + +#endif diff --git a/common/Database.hpp b/common/Database.hpp new file mode 100644 index 00000000..37a89af0 --- /dev/null +++ b/common/Database.hpp @@ -0,0 +1,173 @@ +#ifndef included_Database_hpp +#define included_Database_hpp + +#include "common/Database.h" +#include "common/Utilities.h" + +#include + + +/******************************************************************** + * Basic classes for primative data types * + ********************************************************************/ +class EmptyKeyData : public KeyData +{ +public: + EmptyKeyData() {} + virtual ~EmptyKeyData() {} + virtual std::shared_ptr clone() const override + { + return std::make_shared(); + } + virtual void print( std::ostream& os, const std::string& = "" ) const override + { + os << std::endl; + } + virtual std::string type() const override { return ""; } +}; +class KeyDataDouble : public KeyData +{ +public: + KeyDataDouble() {} + explicit KeyDataDouble( const std::vector& data, const Units& unit ) + : d_data( data ), d_unit( unit ) + { + } + virtual ~KeyDataDouble() {} + virtual std::shared_ptr clone() const override + { + return std::make_shared( d_data, d_unit ); + } + virtual void print( std::ostream& os, const std::string& indent = "" ) const override; + virtual std::string type() const override { return "double"; } + + static std::tuple read( const std::string& ); + static double readValue( const std::string& ); + +public: + std::vector d_data; + Units d_unit; +}; +class KeyDataBool : public KeyData +{ +public: + KeyDataBool() {} + explicit KeyDataBool( const std::vector& data ) : d_data( data ) {} + virtual ~KeyDataBool() {} + virtual std::shared_ptr clone() const override + { + return std::make_shared( d_data ); + } + virtual void print( std::ostream& os, const std::string& indent = "" ) const override + { + os << indent; + for ( size_t i = 0; i < d_data.size(); i++ ) { + if ( i > 0 ) { + os << ", "; + } + if ( d_data[i] ) { + os << "true"; + } else { + os << "false"; + } + } + os << std::endl; + } + virtual std::string type() const override { return "bool"; } + std::vector d_data; +}; +class KeyDataString : public KeyData +{ +public: + KeyDataString() {} + explicit KeyDataString( const std::vector& data ) : d_data( data ) {} + virtual ~KeyDataString() {} + virtual std::shared_ptr clone() const override + { + return std::make_shared( d_data ); + } + virtual void print( std::ostream& os, const std::string& indent = "" ) const override + { + os << indent; + for ( size_t i = 0; i < d_data.size(); i++ ) { + if ( i > 0 ) { + os << ", "; + } + os << '"' << d_data[i] << '"'; + } + os << std::endl; + } + virtual std::string type() const override { return "string"; } + std::vector d_data; +}; + + +/******************************************************************** + * Get a vector * + ********************************************************************/ +template +inline std::vector Database::getVector( + const std::string& key, const std::string& unit ) const +{ + return getVector( key, Units( unit ) ); +} +template +inline void Database::putVector( + const std::string& key, const std::vector& data, const std::string& unit ) +{ + putVector( key, data, Units( unit ) ); +} + + +/******************************************************************** + * Get a scalar * + ********************************************************************/ +template +inline TYPE Database::getScalar( const std::string& key, const Units& unit ) const +{ + const std::vector& data = getVector( key, unit ); + if ( data.size() != 1 ) { + char msg[1000]; + sprintf( msg, "Variable %s is not a scalar", key.c_str() ); + ERROR( msg ); + } + return data[0]; +} +template +inline TYPE Database::getWithDefault( + const std::string& key, const TYPE& value, const Units& unit ) const +{ + if ( !keyExists( key ) ) + return value; + return getScalar( key, unit ); +} +template +inline void Database::putScalar( const std::string& key, const TYPE& data, const Units& unit ) +{ + putVector( key, std::vector( 1, data ), unit ); +} +template +inline TYPE Database::getScalar( const std::string& key, const std::string& unit ) const +{ + return getScalar( key, Units( unit ) ); +} +template +inline TYPE Database::getWithDefault( + const std::string& key, const TYPE& value, const std::string& unit ) const +{ + return getWithDefault( key, value, Units( unit ) ); +} +template +inline void Database::putScalar( const std::string& key, const TYPE& data, const std::string& unit ) +{ + putScalar( key, data, Units( unit ) ); +} +template +inline void putVector( + const std::string& key, const std::vector& data, const std::string& unit ) +{ + putVector( key, data, Units( unit ) ); +} + + +#endif diff --git a/common/Domain.cpp b/common/Domain.cpp index 4c59e2fa..20d8f1fd 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -15,92 +15,131 @@ #include "common/MPI_Helpers.h" #include "common/Communication.h" -using namespace std; -// Reading the domain information file -void read_domain( int rank, int nprocs, MPI_Comm comm, - int& nprocx, int& nprocy, int& nprocz, int& nx, int& ny, int& nz, - int& nspheres, double& Lx, double& Ly, double& Lz ) + + +// Inline function to read line without a return argument +static inline void fgetl( char * str, int num, FILE * stream ) { - if (rank==0){ - ifstream domain("Domain.in"); - domain >> nprocx; - domain >> nprocy; - domain >> nprocz; - domain >> nx; - domain >> ny; - domain >> nz; - domain >> nspheres; - domain >> Lx; - domain >> Ly; - domain >> Lz; - - } - MPI_Barrier(comm); - // Computational domain - //................................................. - MPI_Bcast(&nx,1,MPI_INT,0,comm); - MPI_Bcast(&ny,1,MPI_INT,0,comm); - MPI_Bcast(&nz,1,MPI_INT,0,comm); - MPI_Bcast(&nprocx,1,MPI_INT,0,comm); - MPI_Bcast(&nprocy,1,MPI_INT,0,comm); - MPI_Bcast(&nprocz,1,MPI_INT,0,comm); - MPI_Bcast(&nspheres,1,MPI_INT,0,comm); - MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm); - MPI_Barrier(comm); + char* ptr = fgets( str, num, stream ); + if ( 0 ) {char *temp = (char *)&ptr; temp++;} } + /******************************************************** - * Constructor/Destructor * + * Constructors/Destructor * ********************************************************/ -Domain::Domain(int nx, int ny, int nz, int rnk, int npx, int npy, int npz, - double lx, double ly, double lz, int BC): - Nx(0), Ny(0), Nz(0), iproc(0), jproc(0), nprocx(0), nprocy(0), nprocz(0), - Lx(0), Ly(0), Lz(0), Volume(0), rank(0), BoundaryCondition(0), - Group(MPI_GROUP_NULL), Comm(MPI_COMM_NULL), - rank_x(0), rank_y(0), rank_z(0), rank_X(0), rank_Y(0), rank_Z(0), - rank_xy(0), rank_XY(0), rank_xY(0), rank_Xy(0), - rank_xz(0), rank_XZ(0), rank_xZ(0), rank_Xz(0), - rank_yz(0), rank_YZ(0), rank_yZ(0), rank_Yz(0), - sendCount_x(0), sendCount_y(0), sendCount_z(0), sendCount_X(0), sendCount_Y(0), sendCount_Z(0), - sendCount_xy(0), sendCount_yz(0), sendCount_xz(0), sendCount_Xy(0), sendCount_Yz(0), sendCount_xZ(0), - sendCount_xY(0), sendCount_yZ(0), sendCount_Xz(0), sendCount_XY(0), sendCount_YZ(0), sendCount_XZ(0), - sendList_x(NULL), sendList_y(NULL), sendList_z(NULL), sendList_X(NULL), sendList_Y(NULL), sendList_Z(NULL), - sendList_xy(NULL), sendList_yz(NULL), sendList_xz(NULL), sendList_Xy(NULL), sendList_Yz(NULL), sendList_xZ(NULL), - sendList_xY(NULL), sendList_yZ(NULL), sendList_Xz(NULL), sendList_XY(NULL), sendList_YZ(NULL), sendList_XZ(NULL), - sendBuf_x(NULL), sendBuf_y(NULL), sendBuf_z(NULL), sendBuf_X(NULL), sendBuf_Y(NULL), sendBuf_Z(NULL), - sendBuf_xy(NULL), sendBuf_yz(NULL), sendBuf_xz(NULL), sendBuf_Xy(NULL), sendBuf_Yz(NULL), sendBuf_xZ(NULL), - sendBuf_xY(NULL), sendBuf_yZ(NULL), sendBuf_Xz(NULL), sendBuf_XY(NULL), sendBuf_YZ(NULL), sendBuf_XZ(NULL), - recvCount_x(0), recvCount_y(0), recvCount_z(0), recvCount_X(0), recvCount_Y(0), recvCount_Z(0), - recvCount_xy(0), recvCount_yz(0), recvCount_xz(0), recvCount_Xy(0), recvCount_Yz(0), recvCount_xZ(0), - recvCount_xY(0), recvCount_yZ(0), recvCount_Xz(0), recvCount_XY(0), recvCount_YZ(0), recvCount_XZ(0), - recvList_x(NULL), recvList_y(NULL), recvList_z(NULL), recvList_X(NULL), recvList_Y(NULL), recvList_Z(NULL), - recvList_xy(NULL), recvList_yz(NULL), recvList_xz(NULL), recvList_Xy(NULL), recvList_Yz(NULL), recvList_xZ(NULL), - recvList_xY(NULL), recvList_yZ(NULL), recvList_Xz(NULL), recvList_XY(NULL), recvList_YZ(NULL), recvList_XZ(NULL), - recvBuf_x(NULL), recvBuf_y(NULL), recvBuf_z(NULL), recvBuf_X(NULL), recvBuf_Y(NULL), recvBuf_Z(NULL), - recvBuf_xy(NULL), recvBuf_yz(NULL), recvBuf_xz(NULL), recvBuf_Xy(NULL), recvBuf_Yz(NULL), recvBuf_xZ(NULL), - recvBuf_xY(NULL), recvBuf_yZ(NULL), recvBuf_Xz(NULL), recvBuf_XY(NULL), recvBuf_YZ(NULL), recvBuf_XZ(NULL), - sendData_x(NULL), sendData_y(NULL), sendData_z(NULL), sendData_X(NULL), sendData_Y(NULL), sendData_Z(NULL), - sendData_xy(NULL), sendData_yz(NULL), sendData_xz(NULL), sendData_Xy(NULL), sendData_Yz(NULL), sendData_xZ(NULL), - sendData_xY(NULL), sendData_yZ(NULL), sendData_Xz(NULL), sendData_XY(NULL), sendData_YZ(NULL), sendData_XZ(NULL), - recvData_x(NULL), recvData_y(NULL), recvData_z(NULL), recvData_X(NULL), recvData_Y(NULL), recvData_Z(NULL), - recvData_xy(NULL), recvData_yz(NULL), recvData_xz(NULL), recvData_Xy(NULL), recvData_Yz(NULL), recvData_xZ(NULL), - recvData_xY(NULL), recvData_yZ(NULL), recvData_Xz(NULL), recvData_XY(NULL), recvData_YZ(NULL), recvData_XZ(NULL), - id(NULL) +Domain::Domain( int nx, int ny, int nz, int rnk, int npx, int npy, int npz, + double lx, double ly, double lz, int BC): + Nx(0), Ny(0), Nz(0), iproc(0), jproc(0), nprocx(0), nprocy(0), nprocz(0), + Lx(0), Ly(0), Lz(0), Volume(0), rank(0), BoundaryCondition(0), + Group(MPI_GROUP_NULL), Comm(MPI_COMM_NULL), + rank_x(0), rank_y(0), rank_z(0), rank_X(0), rank_Y(0), rank_Z(0), + rank_xy(0), rank_XY(0), rank_xY(0), rank_Xy(0), + rank_xz(0), rank_XZ(0), rank_xZ(0), rank_Xz(0), + rank_yz(0), rank_YZ(0), rank_yZ(0), rank_Yz(0), + sendCount_x(0), sendCount_y(0), sendCount_z(0), sendCount_X(0), sendCount_Y(0), sendCount_Z(0), + sendCount_xy(0), sendCount_yz(0), sendCount_xz(0), sendCount_Xy(0), sendCount_Yz(0), sendCount_xZ(0), + sendCount_xY(0), sendCount_yZ(0), sendCount_Xz(0), sendCount_XY(0), sendCount_YZ(0), sendCount_XZ(0), + sendList_x(NULL), sendList_y(NULL), sendList_z(NULL), sendList_X(NULL), sendList_Y(NULL), sendList_Z(NULL), + sendList_xy(NULL), sendList_yz(NULL), sendList_xz(NULL), sendList_Xy(NULL), sendList_Yz(NULL), sendList_xZ(NULL), + sendList_xY(NULL), sendList_yZ(NULL), sendList_Xz(NULL), sendList_XY(NULL), sendList_YZ(NULL), sendList_XZ(NULL), + sendBuf_x(NULL), sendBuf_y(NULL), sendBuf_z(NULL), sendBuf_X(NULL), sendBuf_Y(NULL), sendBuf_Z(NULL), + sendBuf_xy(NULL), sendBuf_yz(NULL), sendBuf_xz(NULL), sendBuf_Xy(NULL), sendBuf_Yz(NULL), sendBuf_xZ(NULL), + sendBuf_xY(NULL), sendBuf_yZ(NULL), sendBuf_Xz(NULL), sendBuf_XY(NULL), sendBuf_YZ(NULL), sendBuf_XZ(NULL), + recvCount_x(0), recvCount_y(0), recvCount_z(0), recvCount_X(0), recvCount_Y(0), recvCount_Z(0), + recvCount_xy(0), recvCount_yz(0), recvCount_xz(0), recvCount_Xy(0), recvCount_Yz(0), recvCount_xZ(0), + recvCount_xY(0), recvCount_yZ(0), recvCount_Xz(0), recvCount_XY(0), recvCount_YZ(0), recvCount_XZ(0), + recvList_x(NULL), recvList_y(NULL), recvList_z(NULL), recvList_X(NULL), recvList_Y(NULL), recvList_Z(NULL), + recvList_xy(NULL), recvList_yz(NULL), recvList_xz(NULL), recvList_Xy(NULL), recvList_Yz(NULL), recvList_xZ(NULL), + recvList_xY(NULL), recvList_yZ(NULL), recvList_Xz(NULL), recvList_XY(NULL), recvList_YZ(NULL), recvList_XZ(NULL), + recvBuf_x(NULL), recvBuf_y(NULL), recvBuf_z(NULL), recvBuf_X(NULL), recvBuf_Y(NULL), recvBuf_Z(NULL), + recvBuf_xy(NULL), recvBuf_yz(NULL), recvBuf_xz(NULL), recvBuf_Xy(NULL), recvBuf_Yz(NULL), recvBuf_xZ(NULL), + recvBuf_xY(NULL), recvBuf_yZ(NULL), recvBuf_Xz(NULL), recvBuf_XY(NULL), recvBuf_YZ(NULL), recvBuf_XZ(NULL), + sendData_x(NULL), sendData_y(NULL), sendData_z(NULL), sendData_X(NULL), sendData_Y(NULL), sendData_Z(NULL), + sendData_xy(NULL), sendData_yz(NULL), sendData_xz(NULL), sendData_Xy(NULL), sendData_Yz(NULL), sendData_xZ(NULL), + sendData_xY(NULL), sendData_yZ(NULL), sendData_Xz(NULL), sendData_XY(NULL), sendData_YZ(NULL), sendData_XZ(NULL), + recvData_x(NULL), recvData_y(NULL), recvData_z(NULL), recvData_X(NULL), recvData_Y(NULL), recvData_Z(NULL), + recvData_xy(NULL), recvData_yz(NULL), recvData_xz(NULL), recvData_Xy(NULL), recvData_Yz(NULL), recvData_xZ(NULL), + recvData_xY(NULL), recvData_yZ(NULL), recvData_Xz(NULL), recvData_XY(NULL), recvData_YZ(NULL), recvData_XZ(NULL), + id(NULL) { - Volume = nx*ny*nx*npx*npy*npz*1.0; - Nx = nx+2; Ny = ny+2; Nz = nz+2; - Lx = lx, Ly = ly, Lz = lz; - rank = rnk; - nprocx=npx; nprocy=npy; nprocz=npz; + auto db = std::make_shared( ); + db->putScalar( "BC", BC ); + db->putVector( "nproc", { npx, npx, npx } ); + db->putVector( "n", { nx, ny, nz } ); + db->putScalar( "nspheres", 0 ); + db->putVector( "L", { lx, ly, lz } ); + initialize( db ); +} +Domain::Domain( std::shared_ptr db ): + Nx(0), Ny(0), Nz(0), iproc(0), jproc(0), nprocx(0), nprocy(0), nprocz(0), + Lx(0), Ly(0), Lz(0), Volume(0), rank(0), BoundaryCondition(0), + Group(MPI_GROUP_NULL), Comm(MPI_COMM_NULL), + rank_x(0), rank_y(0), rank_z(0), rank_X(0), rank_Y(0), rank_Z(0), + rank_xy(0), rank_XY(0), rank_xY(0), rank_Xy(0), + rank_xz(0), rank_XZ(0), rank_xZ(0), rank_Xz(0), + rank_yz(0), rank_YZ(0), rank_yZ(0), rank_Yz(0), + sendCount_x(0), sendCount_y(0), sendCount_z(0), sendCount_X(0), sendCount_Y(0), sendCount_Z(0), + sendCount_xy(0), sendCount_yz(0), sendCount_xz(0), sendCount_Xy(0), sendCount_Yz(0), sendCount_xZ(0), + sendCount_xY(0), sendCount_yZ(0), sendCount_Xz(0), sendCount_XY(0), sendCount_YZ(0), sendCount_XZ(0), + sendList_x(NULL), sendList_y(NULL), sendList_z(NULL), sendList_X(NULL), sendList_Y(NULL), sendList_Z(NULL), + sendList_xy(NULL), sendList_yz(NULL), sendList_xz(NULL), sendList_Xy(NULL), sendList_Yz(NULL), sendList_xZ(NULL), + sendList_xY(NULL), sendList_yZ(NULL), sendList_Xz(NULL), sendList_XY(NULL), sendList_YZ(NULL), sendList_XZ(NULL), + sendBuf_x(NULL), sendBuf_y(NULL), sendBuf_z(NULL), sendBuf_X(NULL), sendBuf_Y(NULL), sendBuf_Z(NULL), + sendBuf_xy(NULL), sendBuf_yz(NULL), sendBuf_xz(NULL), sendBuf_Xy(NULL), sendBuf_Yz(NULL), sendBuf_xZ(NULL), + sendBuf_xY(NULL), sendBuf_yZ(NULL), sendBuf_Xz(NULL), sendBuf_XY(NULL), sendBuf_YZ(NULL), sendBuf_XZ(NULL), + recvCount_x(0), recvCount_y(0), recvCount_z(0), recvCount_X(0), recvCount_Y(0), recvCount_Z(0), + recvCount_xy(0), recvCount_yz(0), recvCount_xz(0), recvCount_Xy(0), recvCount_Yz(0), recvCount_xZ(0), + recvCount_xY(0), recvCount_yZ(0), recvCount_Xz(0), recvCount_XY(0), recvCount_YZ(0), recvCount_XZ(0), + recvList_x(NULL), recvList_y(NULL), recvList_z(NULL), recvList_X(NULL), recvList_Y(NULL), recvList_Z(NULL), + recvList_xy(NULL), recvList_yz(NULL), recvList_xz(NULL), recvList_Xy(NULL), recvList_Yz(NULL), recvList_xZ(NULL), + recvList_xY(NULL), recvList_yZ(NULL), recvList_Xz(NULL), recvList_XY(NULL), recvList_YZ(NULL), recvList_XZ(NULL), + recvBuf_x(NULL), recvBuf_y(NULL), recvBuf_z(NULL), recvBuf_X(NULL), recvBuf_Y(NULL), recvBuf_Z(NULL), + recvBuf_xy(NULL), recvBuf_yz(NULL), recvBuf_xz(NULL), recvBuf_Xy(NULL), recvBuf_Yz(NULL), recvBuf_xZ(NULL), + recvBuf_xY(NULL), recvBuf_yZ(NULL), recvBuf_Xz(NULL), recvBuf_XY(NULL), recvBuf_YZ(NULL), recvBuf_XZ(NULL), + sendData_x(NULL), sendData_y(NULL), sendData_z(NULL), sendData_X(NULL), sendData_Y(NULL), sendData_Z(NULL), + sendData_xy(NULL), sendData_yz(NULL), sendData_xz(NULL), sendData_Xy(NULL), sendData_Yz(NULL), sendData_xZ(NULL), + sendData_xY(NULL), sendData_yZ(NULL), sendData_Xz(NULL), sendData_XY(NULL), sendData_YZ(NULL), sendData_XZ(NULL), + recvData_x(NULL), recvData_y(NULL), recvData_z(NULL), recvData_X(NULL), recvData_Y(NULL), recvData_Z(NULL), + recvData_xy(NULL), recvData_yz(NULL), recvData_xz(NULL), recvData_Xy(NULL), recvData_Yz(NULL), recvData_xZ(NULL), + recvData_xY(NULL), recvData_yZ(NULL), recvData_Xz(NULL), recvData_XY(NULL), recvData_YZ(NULL), recvData_XZ(NULL), + id(NULL) +{ + initialize( db ); +} +void Domain::initialize( std::shared_ptr db ) +{ + d_db = db; + auto nproc = d_db->getVector("nproc"); + auto n = d_db->getVector("n"); + auto L = d_db->getVector("L"); + ASSERT( n.size() == 3u ); + ASSERT( L.size() == 3u ); + ASSERT( nproc.size() == 3u ); + int nx = n[0]; + int ny = n[1]; + int nz = n[2]; + Lx = L[0]; + Ly = L[1]; + Lz = L[2]; + Nx = nx+2; + Ny = ny+2; + Nz = nz+2; + nprocx = nproc[0]; + nprocy = nproc[1]; + nprocz = nproc[2]; N = Nx*Ny*Nz; + Volume = nx*ny*nx*nprocx*nprocy*nprocz*1.0; id = new char[N]; memset(id,0,N); - BoundaryCondition = BC; + BoundaryCondition = d_db->getScalar("BC"); + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); rank_info=RankInfoStruct(rank,nprocx,nprocy,nprocz); + int nprocs; + MPI_Comm_size( MPI_COMM_WORLD, &nprocs ); + INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!"); } Domain::~Domain() { @@ -1240,3 +1279,447 @@ void Domain::CommunicateMeshHalo(DoubleArray &Mesh) UnpackMeshData(recvList_YZ, recvCount_YZ ,recvData_YZ, MeshData); } + + + + +double SSO(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){ + /* + * This routine converts the data in the Distance array to a signed distance + * by solving the equation df/dt = sign(1-|grad f|), where Distance provides + * the values of f on the mesh associated with domain Dm + * It has been tested with segmented data initialized to values [-1,1] + * and will converge toward the signed distance to the surface bounding the associated phases + */ + + int Q=26; + int q,i,j,k; + double dt=0.1; + int in,jn,kn; + double Dqx,Dqy,Dqz,Dx,Dy,Dz,W; + double nx,ny,nz,Cqx,Cqy,Cqz,sign,norm; + double TotalVariation=0.0; + + const static int D3Q27[26][3]={{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}, + {1,1,0},{-1,-1,0},{1,-1,0},{-1,1,0},{1,0,1},{-1,0,-1},{1,0,-1},{-1,0,1}, + {0,1,1},{0,-1,-1},{0,1,-1},{0,-1,1},{1,1,1},{-1,-1,-1},{1,1,-1},{-1,-1,1}, + {-1,1,-1},{1,-1,1},{1,-1,-1},{-1,1,1}}; + + double weights[26]; + // Compute the weights from the finite differences + for (q=0; q fillData(Dm.Comm, Dm.rank_info,xdim,ydim,zdim,1,1,1,0,1); + + int count = 0; + while (count < timesteps){ + + // Communicate the halo of values + fillData.fill(Distance); + + TotalVariation=0.0; + // Execute the next timestep + for (k=1;k 0.0 ){ + for (q=0; q<26; q++){ + Cqx = 1.0*D3Q27[q][0]; + Cqy = 1.0*D3Q27[q][1]; + Cqz = 1.0*D3Q27[q][2]; + // get the associated neighbor + in = i + D3Q27[q][0]; + jn = j + D3Q27[q][1]; + kn = k + D3Q27[q][2]; + + // make sure the neighbor is in the domain (periodic BC) + /* if (in < 0 ) in +=Nx; + * don't need this in parallel since MPI handles the halos + if (jn < 0 ) jn +=Ny; + if (kn < 0 ) kn +=Nz; + if (!(in < Nx) ) in -=Nx; + if (!(jn < Ny) ) jn -=Ny; + if (!(kn < Nz) ) kn -=Nz; + // symmetric boundary + if (in < 0 ) in = i; + if (jn < 0 ) jn = j; + if (kn < 0 ) kn = k; + if (!(in < Nx) ) in = i; + if (!(jn < Ny) ) jn = k; + if (!(kn < Nz) ) kn = k; + */ + + // Compute the gradient using upwind finite differences + Dqx = weights[q]*(Distance(i,j,k) - Distance(in,jn,kn))*Cqx; + Dqy = weights[q]*(Distance(i,j,k) - Distance(in,jn,kn))*Cqy; + Dqz = weights[q]*(Distance(i,j,k) - Distance(in,jn,kn))*Cqz; + + // Only include upwind derivatives + if (sign*(nx*Cqx + ny*Cqy + nz*Cqz) < 0.0 ){ + + Dx += Dqx; + Dy += Dqy; + Dz += Dqz; + W += weights[q]; + } + } + // Normalize by the weight to get the approximation to the gradient + if (fabs(W) > 0.0){ + Dx /= W; + Dy /= W; + Dz /= W; + } + norm = sqrt(Dx*Dx+Dy*Dy+Dz*Dz); + } + else{ + norm = 0.0; + } + Distance(i,j,k) += dt*sign*(1.0 - norm); + TotalVariation += dt*sign*(1.0 - norm); + // Disallow any change in phase + // if (Distance(i,j,k)*2.0*(ID[n]-1.0) < 0) Distance(i,j,k) = -Distance(i,j,k); + } + } + } + TotalVariation /= (Dm.Nx-2)*(Dm.Ny-2)*(Dm.Nz-2); + count++; + } + + return TotalVariation; +} + + +void ReadSpherePacking(int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad) +{ + // Read in the full sphere pack + //...... READ IN THE SPHERES................................... + cout << "Reading the packing file..." << endl; + FILE *fid = fopen("pack.out","rb"); + INSIST(fid!=NULL,"Error opening pack.out"); + //.........Trash the header lines.......... + char line[100]; + fgetl(line, 100, fid); + fgetl(line, 100, fid); + fgetl(line, 100, fid); + fgetl(line, 100, fid); + fgetl(line, 100, fid); + //........read the spheres.................. + // We will read until a blank like or end-of-file is reached + int count = 0; + while ( !feof(fid) && fgets(line,100,fid)!=NULL ) { + char* line2 = line; + List_cx[count] = strtod(line2,&line2); + List_cy[count] = strtod(line2,&line2); + List_cz[count] = strtod(line2,&line2); + List_rad[count] = strtod(line2,&line2); + count++; + } + cout << "Number of spheres extracted is: " << count << endl; + INSIST( count==nspheres, "Specified number of spheres is probably incorrect!" ); + // ............................................................. +} + +void AssignLocalSolidID(char *ID, int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad, + double Lx, double Ly, double Lz, int Nx, int Ny, int Nz, + int iproc, int jproc, int kproc, int nprocx, int nprocy, int nprocz) +{ + // Use sphere lists to determine which nodes are in porespace + // Write out binary file for nodes + char value; + int N = Nx*Ny*Nz; // Domain size, including the halo + double hx,hy,hz; + double x,y,z; + double cx,cy,cz,r; + int imin,imax,jmin,jmax,kmin,kmax; + int p,i,j,k,n; + //............................................ + double min_x,min_y,min_z; +// double max_x,max_y,max_z; + //............................................ + // Lattice spacing for the entire domain + // It should generally be true that hx=hy=hz + // Otherwise, you will end up with ellipsoids + hx = Lx/(Nx*nprocx-1); + hy = Ly/(Ny*nprocy-1); + hz = Lz/(Nz*nprocz-1); + //............................................ + // Get maximum and minimum for this domain + // Halo is included ! + min_x = double(iproc*Nx-1)*hx; + min_y = double(jproc*Ny-1)*hy; + min_z = double(kproc*Nz-1)*hz; +// max_x = ((iproc+1)*Nx+1)*hx; +// max_y = ((jproc+1)*Ny+1)*hy; +// max_z = ((kproc+1)*Nz+1)*hz; + //............................................ + + //............................................ + // Pre-initialize local ID + for (n=0;nNx) imin = Nx; + if (imax<0) imax = 0; + if (imax>Nx) imax = Nx; + if (jmin<0) jmin = 0; + if (jmin>Ny) jmin = Ny; + if (jmax<0) jmax = 0; + if (jmax>Ny) jmax = Ny; + if (kmin<0) kmin = 0; + if (kmin>Nz) kmin = Nz; + if (kmax<0) kmax = 0; + if (kmax>Nz) kmax = Nz; + // Loop over the domain for this sphere (may be null) + for (i=imin;iNx) imin = Nx; + if (imax<0) imax = 0; + if (imax>Nx) imax = Nx; + if (jmin<0) jmin = 0; + if (jmin>Ny) jmin = Ny; + if (jmax<0) jmax = 0; + if (jmax>Ny) jmax = Ny; + if (kmin<0) kmin = 0; + if (kmin>Nz) kmin = Nz; + if (kmax<0) kmax = 0; + if (kmax>Nz) kmax = Nz; + // Loop over the domain for this sphere (may be null) + for (i=imin;i #include #include #include #include #include -#include // std::exception +#include #include #include "common/Array.h" #include "common/Utilities.h" #include "common/MPI_Helpers.h" #include "common/Communication.h" +#include "common/Database.h" -using namespace std; //! Read the domain information file -void read_domain( int rank, int nprocs, MPI_Comm comm, - int& nprocx, int& nprocy, int& nprocz, int& nx, int& ny, int& nz, - int& nspheres, double& Lx, double& Ly, double& Lz ); +std::shared_ptr read_domain( ); //! Class to hold domain info -struct Domain{ - // Default constructor - Domain(int nx, int ny, int nz, int rnk, int npx, int npy, int npz, - double lx, double ly, double lz, int BC); +class Domain{ +public: + //! Default constructor + Domain( std::shared_ptr db ); - // Destructor + //! Obsolete constructor + Domain( int nx, int ny, int nz, int rnk, int npx, int npy, int npz, + double lx, double ly, double lz, int BC); + + //! Empty constructor + Domain() = delete; + + //! Copy constructor + Domain( const Domain& ) = delete; + + //! Assignment operator + Domain& operator=( const Domain& ) = delete; + + //! Destructor ~Domain(); + //! Get the database + inline std::shared_ptr getDatabase() const { return d_db; } + +private: + + void initialize( std::shared_ptr db ); + + std::shared_ptr d_db; + +public: + // Basic domain information int Nx,Ny,Nz,N; int iproc,jproc,kproc; @@ -114,454 +135,27 @@ private: }; +double SSO(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps); -// Inline function to read line without a return argument -static inline void fgetl( char * str, int num, FILE * stream ) -{ - char* ptr = fgets( str, num, stream ); - if ( 0 ) {char *temp = (char *)&ptr; temp++;} -} +void ReadSpherePacking(int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad); +void AssignLocalSolidID(char *ID, int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad, + double Lx, double Ly, double Lz, int Nx, int Ny, int Nz, + int iproc, int jproc, int kproc, int nprocx, int nprocy, int nprocz); +void SignedDistance(double *Distance, int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad, + double Lx, double Ly, double Lz, int Nx, int Ny, int Nz, + int iproc, int jproc, int kproc, int nprocx, int nprocy, int nprocz); -inline double SSO(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){ - /* - * This routine converts the data in the Distance array to a signed distance - * by solving the equation df/dt = sign(1-|grad f|), where Distance provides - * the values of f on the mesh associated with domain Dm - * It has been tested with segmented data initialized to values [-1,1] - * and will converge toward the signed distance to the surface bounding the associated phases - */ +void WriteLocalSolidID(char *FILENAME, char *ID, int N); - int Q=26; - int q,i,j,k; - double dt=0.1; - int in,jn,kn; - double Dqx,Dqy,Dqz,Dx,Dy,Dz,W; - double nx,ny,nz,Cqx,Cqy,Cqz,sign,norm; - double TotalVariation=0.0; +void WriteLocalSolidDistance(char *FILENAME, double *Distance, int N); - const static int D3Q27[26][3]={{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}, - {1,1,0},{-1,-1,0},{1,-1,0},{-1,1,0},{1,0,1},{-1,0,-1},{1,0,-1},{-1,0,1}, - {0,1,1},{0,-1,-1},{0,1,-1},{0,-1,1},{1,1,1},{-1,-1,-1},{1,1,-1},{-1,-1,1}, - {-1,1,-1},{1,-1,1},{1,-1,-1},{-1,1,1}}; +void WriteCheckpoint(const char *FILENAME, const double *cDen, const double *cfq, int Np); - double weights[26]; - // Compute the weights from the finite differences - for (q=0; q fillData(Dm.Comm, Dm.rank_info,xdim,ydim,zdim,1,1,1,0,1); - - int count = 0; - while (count < timesteps){ - - // Communicate the halo of values - fillData.fill(Distance); - - TotalVariation=0.0; - // Execute the next timestep - for (k=1;k 0.0 ){ - for (q=0; q<26; q++){ - Cqx = 1.0*D3Q27[q][0]; - Cqy = 1.0*D3Q27[q][1]; - Cqz = 1.0*D3Q27[q][2]; - // get the associated neighbor - in = i + D3Q27[q][0]; - jn = j + D3Q27[q][1]; - kn = k + D3Q27[q][2]; - - // make sure the neighbor is in the domain (periodic BC) - /* if (in < 0 ) in +=Nx; - * don't need this in parallel since MPI handles the halos - if (jn < 0 ) jn +=Ny; - if (kn < 0 ) kn +=Nz; - if (!(in < Nx) ) in -=Nx; - if (!(jn < Ny) ) jn -=Ny; - if (!(kn < Nz) ) kn -=Nz; - // symmetric boundary - if (in < 0 ) in = i; - if (jn < 0 ) jn = j; - if (kn < 0 ) kn = k; - if (!(in < Nx) ) in = i; - if (!(jn < Ny) ) jn = k; - if (!(kn < Nz) ) kn = k; - */ - - // Compute the gradient using upwind finite differences - Dqx = weights[q]*(Distance(i,j,k) - Distance(in,jn,kn))*Cqx; - Dqy = weights[q]*(Distance(i,j,k) - Distance(in,jn,kn))*Cqy; - Dqz = weights[q]*(Distance(i,j,k) - Distance(in,jn,kn))*Cqz; - - // Only include upwind derivatives - if (sign*(nx*Cqx + ny*Cqy + nz*Cqz) < 0.0 ){ - - Dx += Dqx; - Dy += Dqy; - Dz += Dqz; - W += weights[q]; - } - } - // Normalize by the weight to get the approximation to the gradient - if (fabs(W) > 0.0){ - Dx /= W; - Dy /= W; - Dz /= W; - } - norm = sqrt(Dx*Dx+Dy*Dy+Dz*Dz); - } - else{ - norm = 0.0; - } - Distance(i,j,k) += dt*sign*(1.0 - norm); - TotalVariation += dt*sign*(1.0 - norm); - // Disallow any change in phase - // if (Distance(i,j,k)*2.0*(ID[n]-1.0) < 0) Distance(i,j,k) = -Distance(i,j,k); - } - } - } - TotalVariation /= (Dm.Nx-2)*(Dm.Ny-2)*(Dm.Nz-2); - count++; - } - - return TotalVariation; -} - - -inline void ReadSpherePacking(int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad) -{ - // Read in the full sphere pack - //...... READ IN THE SPHERES................................... - cout << "Reading the packing file..." << endl; - FILE *fid = fopen("pack.out","rb"); - INSIST(fid!=NULL,"Error opening pack.out"); - //.........Trash the header lines.......... - char line[100]; - fgetl(line, 100, fid); - fgetl(line, 100, fid); - fgetl(line, 100, fid); - fgetl(line, 100, fid); - fgetl(line, 100, fid); - //........read the spheres.................. - // We will read until a blank like or end-of-file is reached - int count = 0; - while ( !feof(fid) && fgets(line,100,fid)!=NULL ) { - char* line2 = line; - List_cx[count] = strtod(line2,&line2); - List_cy[count] = strtod(line2,&line2); - List_cz[count] = strtod(line2,&line2); - List_rad[count] = strtod(line2,&line2); - count++; - } - cout << "Number of spheres extracted is: " << count << endl; - INSIST( count==nspheres, "Specified number of spheres is probably incorrect!" ); - // ............................................................. -} - -inline void AssignLocalSolidID(char *ID, int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad, - double Lx, double Ly, double Lz, int Nx, int Ny, int Nz, - int iproc, int jproc, int kproc, int nprocx, int nprocy, int nprocz) -{ - // Use sphere lists to determine which nodes are in porespace - // Write out binary file for nodes - char value; - int N = Nx*Ny*Nz; // Domain size, including the halo - double hx,hy,hz; - double x,y,z; - double cx,cy,cz,r; - int imin,imax,jmin,jmax,kmin,kmax; - int p,i,j,k,n; - //............................................ - double min_x,min_y,min_z; -// double max_x,max_y,max_z; - //............................................ - // Lattice spacing for the entire domain - // It should generally be true that hx=hy=hz - // Otherwise, you will end up with ellipsoids - hx = Lx/(Nx*nprocx-1); - hy = Ly/(Ny*nprocy-1); - hz = Lz/(Nz*nprocz-1); - //............................................ - // Get maximum and minimum for this domain - // Halo is included ! - min_x = double(iproc*Nx-1)*hx; - min_y = double(jproc*Ny-1)*hy; - min_z = double(kproc*Nz-1)*hz; -// max_x = ((iproc+1)*Nx+1)*hx; -// max_y = ((jproc+1)*Ny+1)*hy; -// max_z = ((kproc+1)*Nz+1)*hz; - //............................................ - - //............................................ - // Pre-initialize local ID - for (n=0;nNx) imin = Nx; - if (imax<0) imax = 0; - if (imax>Nx) imax = Nx; - if (jmin<0) jmin = 0; - if (jmin>Ny) jmin = Ny; - if (jmax<0) jmax = 0; - if (jmax>Ny) jmax = Ny; - if (kmin<0) kmin = 0; - if (kmin>Nz) kmin = Nz; - if (kmax<0) kmax = 0; - if (kmax>Nz) kmax = Nz; - // Loop over the domain for this sphere (may be null) - for (i=imin;iNx) imin = Nx; - if (imax<0) imax = 0; - if (imax>Nx) imax = Nx; - if (jmin<0) jmin = 0; - if (jmin>Ny) jmin = Ny; - if (jmax<0) jmax = 0; - if (jmax>Ny) jmax = Ny; - if (kmin<0) kmin = 0; - if (kmin>Nz) kmin = Nz; - if (kmax<0) kmax = 0; - if (kmax>Nz) kmax = Nz; - // Loop over the domain for this sphere (may be null) - for (i=imin;i +#include +#include + + +constexpr double Units::d_pow10[22]; +constexpr char Units::d_prefixSymbol[]; + + +/******************************************************************** + * Constructors * + ********************************************************************/ +Units::Units() : d_prefix( UnitPrefix::unknown ), d_unit( UnitValue::unknown ) {} +Units::Units( UnitPrefix p, UnitValue u ) : d_prefix( p ), d_unit( u ) {} +Units::Units( const std::string& unit ) + : d_prefix( UnitPrefix::unknown ), d_unit( UnitValue::unknown ) +{ + // Parse the string to get it into a more friendly format + auto tmp = unit; + tmp.erase( std::remove( tmp.begin(), tmp.end(), ' ' ), tmp.end() ); + // Check if the character '-' is present indicating a seperation between the prefix and unit + size_t index = tmp.find( '-' ); + if ( index != std::string::npos ) { + d_prefix = getUnitPrefix( tmp.substr( 0, index ) ); + d_unit = getUnitValue( tmp.substr( index + 1 ) ); + } else { + if ( tmp.size() <= 1 ) { + d_prefix = UnitPrefix::none; + d_unit = getUnitValue( tmp ); + } else if ( tmp.substr( 0, 2 ) == "da" ) { + d_prefix = UnitPrefix::deca; + d_unit = getUnitValue( tmp.substr( 2 ) ); + } else { + d_prefix = getUnitPrefix( tmp.substr( 0, 1 ) ); + d_unit = getUnitValue( tmp.substr( 1 ) ); + if ( d_prefix == UnitPrefix::unknown || d_unit == UnitValue::unknown ) { + d_prefix = UnitPrefix::none; + d_unit = getUnitValue( tmp ); + } + } + } +} + + +/******************************************************************** + * Get prefix * + ********************************************************************/ +Units::UnitPrefix Units::getUnitPrefix( const std::string& str ) noexcept +{ + Units::UnitPrefix value = UnitPrefix::unknown; + if ( str.empty() ) { + value = UnitPrefix::none; + } else if ( str == "yotta" || str == "Y" ) { + value = UnitPrefix::yotta; + } else if ( str == "zetta" || str == "Z" ) { + value = UnitPrefix::zetta; + } else if ( str == "exa" || str == "E" ) { + value = UnitPrefix::exa; + } else if ( str == "peta" || str == "P" ) { + value = UnitPrefix::peta; + } else if ( str == "tera" || str == "T" ) { + value = UnitPrefix::tera; + } else if ( str == "giga" || str == "G" ) { + value = UnitPrefix::giga; + } else if ( str == "mega" || str == "M" ) { + value = UnitPrefix::mega; + } else if ( str == "kilo" || str == "k" ) { + value = UnitPrefix::kilo; + } else if ( str == "hecto" || str == "h" ) { + value = UnitPrefix::hecto; + } else if ( str == "deca" || str == "da" ) { + value = UnitPrefix::deca; + } else if ( str == "deci" || str == "d" ) { + value = UnitPrefix::deci; + } else if ( str == "centi" || str == "c" ) { + value = UnitPrefix::centi; + } else if ( str == "milli" || str == "m" ) { + value = UnitPrefix::milli; + } else if ( str == "micro" || str == "u" ) { + value = UnitPrefix::micro; + } else if ( str == "nano" || str == "n" ) { + value = UnitPrefix::nano; + } else if ( str == "pico" || str == "p" ) { + value = UnitPrefix::pico; + } else if ( str == "femto" || str == "f" ) { + value = UnitPrefix::femto; + } else if ( str == "atto" || str == "a" ) { + value = UnitPrefix::atto; + } else if ( str == "zepto" || str == "z" ) { + value = UnitPrefix::zepto; + } else if ( str == "yocto" || str == "y" ) { + value = UnitPrefix::yocto; + } + return value; +} + + +/******************************************************************** + * Get unit value * + ********************************************************************/ +Units::UnitValue Units::getUnitValue( const std::string& str ) noexcept +{ + Units::UnitValue value = UnitValue::unknown; + if ( str == "meter" || str == "m" ) { + value = UnitValue::meter; + } else if ( str == "gram" || str == "g" ) { + value = UnitValue::gram; + } else if ( str == "second" || str == "s" ) { + value = UnitValue::second; + } else if ( str == "ampere" || str == "A" ) { + value = UnitValue::ampere; + } else if ( str == "kelvin" || str == "K" ) { + value = UnitValue::kelvin; + } else if ( str == "joule" || str == "J" ) { + value = UnitValue::joule; + } else if ( str == "ergs" || str == "erg" ) { + value = UnitValue::erg; + } else if ( str == "degree" || str == "degrees" ) { + value = UnitValue::degree; + } else if ( str == "radian" || str == "radians" ) { + value = UnitValue::radian; + } + return value; +} + + +/******************************************************************** + * Get unit type * + ********************************************************************/ +Units::UnitType Units::getUnitType( UnitValue u ) noexcept +{ + switch ( u ) { + case UnitValue::meter: + return UnitType::length; + case UnitValue::gram: + return UnitType::mass; + case UnitValue::second: + return UnitType::time; + case UnitValue::ampere: + return UnitType::current; + case UnitValue::kelvin: + return UnitType::temperature; + case UnitValue::joule: + case UnitValue::erg: + return UnitType::energy; + case UnitValue::degree: + case UnitValue::radian: + return UnitType::angle; + default: + return UnitType::unknown; + } +} + + +/******************************************************************** + * Convert to another unit system * + ********************************************************************/ +double Units::convert( const Units& rhs ) const noexcept +{ + if ( this->operator==( rhs ) ) + return 1; + // Convert the prefix + double cp = convert( d_prefix ) / convert( rhs.d_prefix ); + if ( d_unit == rhs.d_unit ) + return cp; // Only need to convert prefix + // Convert the unit + if ( getUnitType( d_unit ) != getUnitType( rhs.d_unit ) ) + return 0; // Invalid conversion + double cu = 0; + if ( d_unit == UnitValue::joule && rhs.d_unit == UnitValue::erg ) + cu = 1e7; + else if ( d_unit == UnitValue::erg && rhs.d_unit == UnitValue::joule ) + cu = 1e-7; + else if ( d_unit == UnitValue::degree && rhs.d_unit == UnitValue::radian ) + cu = 0.017453292519943; + else if ( d_unit == UnitValue::radian && rhs.d_unit == UnitValue::degree ) + cu = 57.295779513082323; + // Return the total conversion + return cp * cu; +} + + +/******************************************************************** + * Write a string for the units * + ********************************************************************/ +std::string Units::str() const +{ + ASSERT( !isNull() ); + return std::string( str( d_prefix ).data() ) + str( d_unit ); +} +std::array Units::str( UnitPrefix p ) noexcept +{ + std::array str; + str[0] = d_prefixSymbol[static_cast( p )]; + str[1] = 0; + str[2] = 0; + if ( p == UnitPrefix::deca ) + str[1] = 'a'; + return str; +} +std::string Units::str( UnitValue u ) +{ + if ( u == UnitValue::meter ) { + return "m"; + } else if ( u == UnitValue::gram ) { + return "g"; + } else if ( u == UnitValue::second ) { + return "s"; + } else if ( u == UnitValue::ampere ) { + return "A"; + } else if ( u == UnitValue::kelvin ) { + return "K"; + } else if ( u == UnitValue::joule ) { + return "J"; + } else if ( u == UnitValue::erg ) { + return "erg"; + } else if ( u == UnitValue::degree ) { + return "degree"; + } else if ( u == UnitValue::radian ) { + return "radian"; + } + return "unknown"; +} diff --git a/common/Units.h b/common/Units.h new file mode 100644 index 00000000..a3496ee1 --- /dev/null +++ b/common/Units.h @@ -0,0 +1,140 @@ +#ifndef included_Units +#define included_Units + +#include +#include +#include +#include +#include +#include + + +//! Unit system class +class Units final +{ +public: + //! Enum to hold prefix + enum class UnitPrefix : int8_t { + yocto = 0, + zepto = 1, + atto = 2, + femto = 3, + pico = 4, + nano = 5, + micro = 6, + milli = 7, + centi = 8, + deci = 9, + none = 10, + deca = 11, + hecto = 12, + kilo = 13, + mega = 14, + giga = 15, + tera = 16, + peta = 17, + exa = 18, + zetta = 19, + yotta = 20, + unknown = 21 + }; + + //! Enum to unit type + enum class UnitType : uint8_t { + length, + mass, + time, + current, + temperature, + energy, + angle, + unknown + }; + + //! Enum to hold unit + enum class UnitValue : uint8_t { + meter, + gram, + second, + ampere, + kelvin, + joule, + erg, + degree, + radian, + unknown + }; + + +public: + //! Constructor + Units(); + + //! Constructor + explicit Units( const std::string& unit ); + + //! Constructor + explicit Units( UnitPrefix, UnitValue ); + + //! Get the prefix + inline UnitPrefix getPrefix() const noexcept { return d_prefix; } + + //! Get the unit + inline UnitValue getUnit() const noexcept { return d_unit; } + + //! Get the unit + inline UnitType getUnitType() const noexcept { return getUnitType( d_unit ); } + + //! Get the unit + static UnitType getUnitType( UnitValue ) noexcept; + + //! Get the prefix from a string + static UnitPrefix getUnitPrefix( const std::string& ) noexcept; + + //! Get the unit value from a string + static UnitValue getUnitValue( const std::string& ) noexcept; + + //! Convert to the given unit system + double convert( const Units& ) const noexcept; + + //! Convert a prefix to a scalar + static inline double convert( UnitPrefix x ) noexcept + { + return d_pow10[static_cast( x )]; + } + + //! Get a string representation of the units + std::string str() const; + + //! Get a string representation for the prefix + static std::array str( UnitPrefix ) noexcept; + + //! Get a string representation for the unit value + static std::string str( UnitValue ); + + //! Operator == + inline bool operator==( const Units& rhs ) const noexcept + { + return d_prefix == rhs.d_prefix && d_unit == rhs.d_unit; + } + + //! Operator != + inline bool operator!=( const Units& rhs ) const noexcept + { + return d_prefix != rhs.d_prefix || d_unit != rhs.d_unit; + } + + //! Check if unit is null + bool isNull() const { return d_prefix == UnitPrefix::unknown || d_unit == UnitValue::unknown; } + +protected: + UnitPrefix d_prefix; + UnitValue d_unit; + +private: + constexpr static double d_pow10[22] = { 1e-24, 1e-21, 1e-18, 1e-15, 1e-12, 1e-9, 1e-6, 1e-3, + 1e-2, 0.1, 1, 10, 100, 1000, 1e6, 1e9, 1e12, 1e15, 1e18, 1e21, 1e24, 0 }; + constexpr static char d_prefixSymbol[] = "yzafpnumcd\0dhkMGTPEZYu"; +}; + +#endif diff --git a/tests/TestColorSquareTube.cpp b/tests/TestColorSquareTube.cpp index 39d16e59..c219e9b4 100644 --- a/tests/TestColorSquareTube.cpp +++ b/tests/TestColorSquareTube.cpp @@ -9,7 +9,35 @@ #include "common/ScaLBL.h" #include "common/MPI_Helpers.h" -using namespace std; + +std::shared_ptr loadInputs( int nprocs ) +{ + auto db = std::make_shared( "Domain.in" ); + const int dim = 50; + db->putScalar( "BC", 0 ); + if ( nprocs == 1 ){ + db->putVector( "nproc", { 1, 1, 1 } ); + db->putVector( "n", { 3, 1, 1 } ); + db->putScalar( "nspheres", 0 ); + db->putVector( "L", { 1, 1, 1 } ); + } else if ( nprocs == 2 ) { + db->putVector( "nproc", { 2, 1, 1 } ); + db->putVector( "n", { dim, dim, dim } ); + db->putScalar( "nspheres", 0 ); + db->putVector( "L", { 1, 1, 1 } ); + } else if ( nprocs == 4 ) { + db->putVector( "nproc", { 2, 2, 1 } ); + db->putVector( "n", { dim, dim, dim } ); + db->putScalar( "nspheres", 0 ); + db->putVector( "L", { 1, 1, 1 } ); + } else if (nprocs==8){ + db->putVector( "nproc", { 2, 2, 2 } ); + db->putVector( "n", { dim, dim, dim } ); + db->putScalar( "nspheres", 0 ); + db->putVector( "L", { 1, 1, 1 } ); + } + return db; +} //*************************************************************************************** @@ -26,11 +54,6 @@ int main(int argc, char **argv) MPI_Comm_size(comm,&nprocs); int check; { - // parallel domain size (# of sub-domains) - int nprocx,nprocy,nprocz; - int iproc,jproc,kproc; - - if (rank == 0){ printf("********************************************************\n"); printf("Running Color Model: TestColor \n"); @@ -43,11 +66,8 @@ int main(int argc, char **argv) int timestepMax, interval; double Fx,Fy,Fz,tol; // Domain variables - double Lx,Ly,Lz; - int nspheres; - int Nx,Ny,Nz; int i,j,k,n; - int dim = 50; + //if (rank == 0) printf("dim=%d\n",dim); int timestep = 1; int timesteps = 100; @@ -78,90 +98,25 @@ int main(int argc, char **argv) flux = 10.f; dout=1.f; - if (rank==0){ - //....................................................................... - // Reading the domain information file - //....................................................................... - ifstream domain("Domain.in"); - if (domain.good()){ - domain >> nprocx; - domain >> nprocy; - domain >> nprocz; - domain >> Nx; - domain >> Ny; - domain >> Nz; - domain >> nspheres; - domain >> Lx; - domain >> Ly; - domain >> Lz; - } - else if (nprocs==1){ - nprocx=nprocy=nprocz=1; - Nx=3; Ny = 1; - Nz = 1; - nspheres=0; - Lx=Ly=Lz=1; - } - else if (nprocs==2){ - nprocx=2; nprocy=1; - nprocz=1; - Nx=Ny=Nz=dim; - Nx = dim; Ny = dim; Nz = dim; - nspheres=0; - Lx=Ly=Lz=1; - } - else if (nprocs==4){ - nprocx=nprocy=2; - nprocz=1; - Nx=Ny=Nz=dim; - nspheres=0; - Lx=Ly=Lz=1; - } - else if (nprocs==8){ - nprocx=nprocy=nprocz=2; - Nx=Ny=Nz=dim; - nspheres=0; - Lx=Ly=Lz=1; - } - //....................................................................... - } - // ************************************************************** - // Broadcast simulation parameters from rank 0 to all other procs - MPI_Barrier(comm); - //................................................. - MPI_Bcast(&Nx,1,MPI_INT,0,comm); - MPI_Bcast(&Ny,1,MPI_INT,0,comm); - MPI_Bcast(&Nz,1,MPI_INT,0,comm); - MPI_Bcast(&nprocx,1,MPI_INT,0,comm); - MPI_Bcast(&nprocy,1,MPI_INT,0,comm); - MPI_Bcast(&nprocz,1,MPI_INT,0,comm); - MPI_Bcast(&nspheres,1,MPI_INT,0,comm); - MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm); - //................................................. - MPI_Barrier(comm); - // ************************************************************** - // ************************************************************** - - if (nprocs != nprocx*nprocy*nprocz){ - printf("nprocx = %i \n",nprocx); - printf("nprocy = %i \n",nprocy); - printf("nprocz = %i \n",nprocz); - INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!"); - } + // Load inputs + auto db = loadInputs( nprocs ); + int Nx = db->getVector( "n" )[0]; + int Ny = db->getVector( "n" )[1]; + int Nz = db->getVector( "n" )[2]; + int nprocx = db->getVector( "nproc" )[0]; + int nprocy = db->getVector( "nproc" )[1]; + int nprocz = db->getVector( "nproc" )[2]; if (rank==0){ printf("********************************************************\n"); printf("Sub-domain size = %i x %i x %i\n",Nx,Ny,Nz); printf("********************************************************\n"); } - MPI_Barrier(comm); double iVol_global = 1.0/Nx/Ny/Nz/nprocx/nprocy/nprocz; - int BoundaryCondition=0; - Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); + + Domain Dm(db); Nx += 2; Ny += 2; @@ -186,9 +141,9 @@ int main(int argc, char **argv) } } - kproc = rank/(nprocx*nprocy); - jproc = (rank-nprocx*nprocy*kproc)/nprocx; - iproc = rank-nprocx*nprocy*kproc-nprocx*jproc; + int kproc = rank/(nprocx*nprocy); + int jproc = (rank-nprocx*nprocy*kproc)/nprocx; + int iproc = rank-nprocx*nprocy*kproc-nprocx*jproc; printf("rank=%i, %i,%i,%i \n",rank,iproc,jproc,kproc); // Initialize a square tube for (k=1;k loadInputs( int nprocs ) +{ + auto db = std::make_shared( "Domain.in" ); + const int dim = 50; + db->putScalar( "BC", 0 ); + if ( nprocs == 1 ){ + db->putVector( "nproc", { 1, 1, 1 } ); + db->putVector( "n", { 3, 1, 1 } ); + db->putScalar( "nspheres", 0 ); + db->putVector( "L", { 1, 1, 1 } ); + } else if ( nprocs == 2 ) { + db->putVector( "nproc", { 2, 1, 1 } ); + db->putVector( "n", { dim, dim, dim } ); + db->putScalar( "nspheres", 0 ); + db->putVector( "L", { 1, 1, 1 } ); + } else if ( nprocs == 4 ) { + db->putVector( "nproc", { 2, 2, 1 } ); + db->putVector( "n", { dim, dim, dim } ); + db->putScalar( "nspheres", 0 ); + db->putVector( "L", { 1, 1, 1 } ); + } else if (nprocs==8){ + db->putVector( "nproc", { 2, 2, 2 } ); + db->putVector( "n", { dim, dim, dim } ); + db->putScalar( "nspheres", 0 ); + db->putVector( "L", { 1, 1, 1 } ); + } + return db; +} + + extern void GlobalFlipScaLBL_D3Q19_Init(double *dist, IntArray Map, int Np, int Nx, int Ny, int Nz, int iproc, int jproc, int kproc, int nprocx, int nprocy, int nprocz) { @@ -176,88 +206,13 @@ int main(int argc, char **argv) int timestepMax, interval; double tau,Fx,Fy,Fz,tol; // Domain variables - double Lx,Ly,Lz; - int nspheres; - int Nx,Ny,Nz; int i,j,k,n; - if (rank==0){ - //....................................................................... - // Reading the domain information file - //....................................................................... - ifstream domain("Domain.in"); - if (domain.good()){ - domain >> nprocx; - domain >> nprocy; - domain >> nprocz; - domain >> Nx; - domain >> Ny; - domain >> Nz; - domain >> nspheres; - domain >> Lx; - domain >> Ly; - domain >> Lz; - } - else if (nprocs==1){ - nprocx=nprocy=nprocz=1; - Nx=Ny=Nz=50; - nspheres=0; - Lx=Ly=Lz=1; - } - else if (nprocs==2){ - nprocx=nprocy=1; - nprocz=2; - Nx=Ny=Nz=50; - nspheres=0; - Lx=Ly=Lz=1; - } - else if (nprocs==4){ - nprocx=nprocy=2; - nprocz=1; - Nx=Ny=Nz=50; - nspheres=0; - Lx=Ly=Lz=1; - } - else if (nprocs==8){ - nprocx=nprocy=nprocz=2; - Nx=Ny=Nz=50; - nspheres=0; - Lx=Ly=Lz=1; - } - //....................................................................... - } - // ************************************************************** - // Broadcast simulation parameters from rank 0 to all other procs - MPI_Barrier(comm); - //................................................. - MPI_Bcast(&Nx,1,MPI_INT,0,comm); - MPI_Bcast(&Ny,1,MPI_INT,0,comm); - MPI_Bcast(&Nz,1,MPI_INT,0,comm); - MPI_Bcast(&nBlocks,1,MPI_UNSIGNED,0,comm); - MPI_Bcast(&nthreads,1,MPI_UNSIGNED,0,comm); - MPI_Bcast(×tepMax,1,MPI_INT,0,comm); - - MPI_Bcast(&Nx,1,MPI_INT,0,comm); - MPI_Bcast(&Ny,1,MPI_INT,0,comm); - MPI_Bcast(&Nz,1,MPI_INT,0,comm); - MPI_Bcast(&nprocx,1,MPI_INT,0,comm); - MPI_Bcast(&nprocy,1,MPI_INT,0,comm); - MPI_Bcast(&nprocz,1,MPI_INT,0,comm); - MPI_Bcast(&nspheres,1,MPI_INT,0,comm); - MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm); - //................................................. - MPI_Barrier(comm); - // ************************************************************** - // ************************************************************** - - if (nprocs != nprocx*nprocy*nprocz){ - printf("nprocx = %i \n",nprocx); - printf("nprocy = %i \n",nprocy); - printf("nprocz = %i \n",nprocz); - INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!"); - } + // Load inputs + auto db = loadInputs( nprocs ); + int Nx = db->getVector( "n" )[0]; + int Ny = db->getVector( "n" )[1]; + int Nz = db->getVector( "n" )[2]; if (rank==0){ printf("********************************************************\n"); @@ -272,8 +227,9 @@ int main(int argc, char **argv) iproc = rank-nprocx*nprocy*kproc-nprocz*jproc; double iVol_global = 1.0/Nx/Ny/Nz/nprocx/nprocy/nprocz; - int BoundaryCondition=0; - Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); + + Domain Dm(db); + InitializeRanks( rank, nprocx, nprocy, nprocz, iproc, jproc, kproc, rank_x, rank_y, rank_z, rank_X, rank_Y, rank_Z, diff --git a/tests/TestForceMoments.cpp b/tests/TestForceMoments.cpp index d07a7213..fe842dc7 100644 --- a/tests/TestForceMoments.cpp +++ b/tests/TestForceMoments.cpp @@ -30,6 +30,22 @@ extern void PrintNeighborList(int * neighborList, int Np, int rank) { } } + +std::shared_ptr loadInputs( int nprocs ) +{ + auto db = std::make_shared( "Domain.in" ); + const int dim = 50; + db->putScalar( "BC", 0 ); + if ( nprocs == 1 ){ + db->putVector( "nproc", { 1, 1, 1 } ); + db->putVector( "n", { 3, 3, 3 } ); + db->putScalar( "nspheres", 0 ); + db->putVector( "L", { 1, 1, 1 } ); + } + return db; +} + + //*************************************************************************************** int main(int argc, char **argv) { @@ -45,7 +61,6 @@ int main(int argc, char **argv) int check; { // parallel domain size (# of sub-domains) - int nprocx,nprocy,nprocz; int iproc,jproc,kproc; @@ -61,9 +76,6 @@ int main(int argc, char **argv) int timestepMax, interval; double tau,Fx,Fy,Fz,tol; // Domain variables - double Lx,Ly,Lz; - int nspheres; - int Nx,Ny,Nz; int i,j,k,n; int dim = 3; if (rank == 0) printf("dim=%d\n",dim); int timestep = 0; @@ -76,58 +88,15 @@ int main(int argc, char **argv) Fx = Fy = 1.0; Fz = 1.0; - if (rank==0){ - //....................................................................... - // Reading the domain information file - //....................................................................... - ifstream domain("Domain.in"); - if (domain.good()){ - printf("domain.good == true \n"); - domain >> nprocx; - domain >> nprocy; - domain >> nprocz; - domain >> Nx; - domain >> Ny; - domain >> Nz; - domain >> nspheres; - domain >> Lx; - domain >> Ly; - domain >> Lz; - } - else if (nprocs==1){ - nprocx=nprocy=nprocz=1; - Nx=3; Ny = 3; - Nz = 3; - nspheres=0; - Lx=Ly=Lz=1; - } - //....................................................................... - } - // ************************************************************** - // Broadcast simulation parameters from rank 0 to all other procs - MPI_Barrier(comm); - //................................................. - MPI_Bcast(&Nx,1,MPI_INT,0,comm); - MPI_Bcast(&Ny,1,MPI_INT,0,comm); - MPI_Bcast(&Nz,1,MPI_INT,0,comm); - MPI_Bcast(&nprocx,1,MPI_INT,0,comm); - MPI_Bcast(&nprocy,1,MPI_INT,0,comm); - MPI_Bcast(&nprocz,1,MPI_INT,0,comm); - MPI_Bcast(&nspheres,1,MPI_INT,0,comm); - MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm); - //................................................. - MPI_Barrier(comm); - // ************************************************************** - // ************************************************************** - if (nprocs != nprocx*nprocy*nprocz){ - printf("nprocx = %i \n",nprocx); - printf("nprocy = %i \n",nprocy); - printf("nprocz = %i \n",nprocz); - INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!"); - } + // Load inputs + auto db = loadInputs( nprocs ); + int Nx = db->getVector( "n" )[0]; + int Ny = db->getVector( "n" )[1]; + int Nz = db->getVector( "n" )[2]; + int nprocx = db->getVector( "nproc" )[0]; + int nprocy = db->getVector( "nproc" )[1]; + int nprocz = db->getVector( "nproc" )[2]; if (rank==0){ printf("********************************************************\n"); @@ -150,9 +119,8 @@ int main(int argc, char **argv) } double iVol_global = 1.0/Nx/Ny/Nz/nprocx/nprocy/nprocz; - int BoundaryCondition=0; - Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); + Domain Dm(db); Nx += 2; Ny += 2; diff --git a/tests/TestInterfaceSpeed.cpp b/tests/TestInterfaceSpeed.cpp index 9ace6822..53cc55ae 100644 --- a/tests/TestInterfaceSpeed.cpp +++ b/tests/TestInterfaceSpeed.cpp @@ -15,6 +15,20 @@ #define SPEED -1 #define PI 3.14159 + +std::shared_ptr loadInputs( int nprocs ) +{ + auto db = std::make_shared( ); + const int dim = 50; + db->putScalar( "BC", 0 ); + db->putVector( "nproc", { 1, 1, 1 } ); + db->putVector( "n", { N, N, N } ); + db->putScalar( "nspheres", 0 ); + db->putVector( "L", { 1, 1, 1 } ); + return db; +} + + int main (int argc, char *argv[]) { // Initialize MPI @@ -24,16 +38,11 @@ int main (int argc, char *argv[]) MPI_Comm_rank(comm,&rank); MPI_Comm_size(comm,&nprocs); - int npx,npy,npz; int i,j,k,n; - int Nx,Ny,Nz; - double Lx,Ly,Lz; - Nx=Ny=Nz=N; - npx=npy=npz=1; - Lx=Ly=Lz=1.0; - int BC=0; // periodic boundary condition + // Load inputs + auto db = loadInputs( nprocs ); - Domain Dm(Nx,Ny,Nz,rank,npx,npy,npz,Lx,Ly,Lz,BC); + Domain Dm(db); for (i=0; i loadInputs( int nprocs ) +{ + const int dim = 12; + std::shared_ptr db; + if ( exists( "Domain.in" ) ) { + db = std::make_shared( "Domain.in" ); + } else if (nprocs==1) { + db->putVector( "nproc", { 1, 1, 1 } ); + db->putVector( "n", { dim, dim, dim } ); + db->putScalar( "nspheres", 0 ); + db->putVector( "L", { 1, 1, 1 } ); + } else if (nprocs==2) { + db->putVector( "nproc", { 1, 1, 1 } ); + db->putVector( "n", { dim/2, dim/2, dim/2 } ); + db->putScalar( "nspheres", 0 ); + db->putVector( "L", { 1, 1, 1 } ); + } + db->putScalar( "BC", 0 ); + return db; +} + //*************************************************************************************** int main(int argc, char **argv) @@ -25,11 +53,6 @@ int main(int argc, char **argv) MPI_Comm_size(comm,&nprocs); int check; { - // parallel domain size (# of sub-domains) - int nprocx,nprocy,nprocz; - int iproc,jproc,kproc; - - if (rank == 0){ printf("********************************************************\n"); printf("Running Unit Test: TestPoiseuille \n"); @@ -42,11 +65,7 @@ int main(int argc, char **argv) int timestepMax, interval; double tau,Fx,Fy,Fz,tol; // Domain variables - double Lx,Ly,Lz; - int nspheres; - int Nx,Ny,Nz; int i,j,k,n; - int dim = 12; if (rank == 0) printf("dim=%d\n",dim); int timestep = 0; @@ -56,68 +75,16 @@ int main(int argc, char **argv) double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); Fx = 0; Fy = 0; Fz = 1e-3; //1.f; // 1e-3; - if (rank==0){ - //....................................................................... - // Reading the domain information file - //....................................................................... - ifstream domain("Domain.in"); - if (domain.good()){ - printf("domain.good == true \n"); - domain >> nprocx; - domain >> nprocy; - domain >> nprocz; - domain >> Nx; - domain >> Ny; - domain >> Nz; - domain >> nspheres; - domain >> Lx; - domain >> Ly; - domain >> Lz; - } - else if (nprocs==1){ - printf("domain.good == false - using predefined parameters \n"); - nprocx=nprocy=nprocz=1; - Nx=dim; Ny = dim; - Nz = dim; - nspheres=0; - Lx=Ly=Lz=1; - } - else if (nprocs==2){ - printf("domain.good == false - using predefined parameters \n"); - nprocx=2; - nprocy=nprocz=1; - Nx=0.5*dim; Ny = 0.5*dim; - Nz = 0.5*dim; - nspheres=0; - Lx=Ly=Lz=1; - } - //....................................................................... - } - // ************************************************************** - // Broadcast simulation parameters from rank 0 to all other procs - MPI_Barrier(comm); - //................................................. - MPI_Bcast(&Nx,1,MPI_INT,0,comm); - MPI_Bcast(&Ny,1,MPI_INT,0,comm); - MPI_Bcast(&Nz,1,MPI_INT,0,comm); - MPI_Bcast(&nprocx,1,MPI_INT,0,comm); - MPI_Bcast(&nprocy,1,MPI_INT,0,comm); - MPI_Bcast(&nprocz,1,MPI_INT,0,comm); - MPI_Bcast(&nspheres,1,MPI_INT,0,comm); - MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm); - //................................................. - MPI_Barrier(comm); - // ************************************************************** - // ************************************************************** - if (nprocs != nprocx*nprocy*nprocz){ - printf("nprocx = %i \n",nprocx); - printf("nprocy = %i \n",nprocy); - printf("nprocz = %i \n",nprocz); - INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!"); - } + // Load inputs + auto db = loadInputs( nprocs ); + int Nx = db->getVector( "n" )[0]; + int Ny = db->getVector( "n" )[1]; + int Nz = db->getVector( "n" )[2]; + int nprocx = db->getVector( "nproc" )[0]; + int nprocy = db->getVector( "nproc" )[1]; + int nprocz = db->getVector( "nproc" )[2]; + if (rank==0){ printf("********************************************************\n"); @@ -126,9 +93,9 @@ int main(int argc, char **argv) } MPI_Barrier(comm); - kproc = rank/(nprocx*nprocy); - jproc = (rank-nprocx*nprocy*kproc)/nprocx; - iproc = rank-nprocx*nprocy*kproc-nprocz*jproc; + int kproc = rank/(nprocx*nprocy); + int jproc = (rank-nprocx*nprocy*kproc)/nprocx; + int iproc = rank-nprocx*nprocy*kproc-nprocz*jproc; if (rank == 0) { printf("i,j,k proc=%d %d %d \n",iproc,jproc,kproc); @@ -140,9 +107,8 @@ int main(int argc, char **argv) } double iVol_global = 1.0/Nx/Ny/Nz/nprocx/nprocy/nprocz; - int BoundaryCondition=0; - Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); + Domain Dm(db); Nx += 2; Ny += 2; diff --git a/tests/TestPressVel.cpp b/tests/TestPressVel.cpp index 233f185d..bdd67aec 100644 --- a/tests/TestPressVel.cpp +++ b/tests/TestPressVel.cpp @@ -9,7 +9,21 @@ #include "common/ScaLBL.h" #include "common/MPI_Helpers.h" -using namespace std; + +std::shared_ptr loadInputs( int nprocs ) +{ + auto db = std::make_shared( "Domain.in" ); + const int dim = 50; + db->putScalar( "BC", 0 ); + if ( nprocs == 1 ){ + db->putVector( "nproc", { 1, 1, 1 } ); + db->putVector( "n", { 4, 4, 4 } ); + db->putScalar( "nspheres", 0 ); + db->putVector( "L", { 1, 1, 1 } ); + } + return db; +} + //*************************************************************************************** int main(int argc, char **argv) @@ -25,11 +39,6 @@ int main(int argc, char **argv) MPI_Comm_size(comm,&nprocs); int check; { - // parallel domain size (# of sub-domains) - int nprocx,nprocy,nprocz; - int iproc,jproc,kproc; - - if (rank == 0){ printf("********************************************************\n"); printf("Running Unit Test: TestPressVel \n"); @@ -42,67 +51,21 @@ int main(int argc, char **argv) int timestepMax, interval; double tau,Fx,Fy,Fz,tol; // Domain variables - double Lx,Ly,Lz; - int nspheres; - int Nx,Ny,Nz; int i,j,k,n; - int dim = 50; if (rank == 0) printf("dim=%d\n",dim); double rlx_setA=1.0; double rlx_setB=1.0; Fx = Fy = 0.f; Fz = 1.0e-4; - if (rank==0){ - //....................................................................... - // Reading the domain information file - //....................................................................... - ifstream domain("Domain.in"); - if (domain.good()){ - domain >> nprocx; - domain >> nprocy; - domain >> nprocz; - domain >> Nx; - domain >> Ny; - domain >> Nz; - domain >> nspheres; - domain >> Lx; - domain >> Ly; - domain >> Lz; - } - else if (nprocs==1){ - nprocx=nprocy=nprocz=1; - Nx=4; Ny = 4; - Nz = 4; - nspheres=0; - Lx=Ly=Lz=1; - } - //....................................................................... - } - // ************************************************************** - // Broadcast simulation parameters from rank 0 to all other procs - MPI_Barrier(comm); - //................................................. - MPI_Bcast(&Nx,1,MPI_INT,0,comm); - MPI_Bcast(&Ny,1,MPI_INT,0,comm); - MPI_Bcast(&Nz,1,MPI_INT,0,comm); - MPI_Bcast(&nprocx,1,MPI_INT,0,comm); - MPI_Bcast(&nprocy,1,MPI_INT,0,comm); - MPI_Bcast(&nprocz,1,MPI_INT,0,comm); - MPI_Bcast(&nspheres,1,MPI_INT,0,comm); - MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm); - //................................................. - MPI_Barrier(comm); - // ************************************************************** - // ************************************************************** - if (nprocs != nprocx*nprocy*nprocz){ - printf("nprocx = %i \n",nprocx); - printf("nprocy = %i \n",nprocy); - printf("nprocz = %i \n",nprocz); - INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!"); - } + // Load inputs + auto db = loadInputs( nprocs ); + int Nx = db->getVector( "n" )[0]; + int Ny = db->getVector( "n" )[1]; + int Nz = db->getVector( "n" )[2]; + int nprocx = db->getVector( "nproc" )[0]; + int nprocy = db->getVector( "nproc" )[1]; + int nprocz = db->getVector( "nproc" )[2]; if (rank==0){ printf("********************************************************\n"); @@ -111,9 +74,9 @@ int main(int argc, char **argv) } MPI_Barrier(comm); - kproc = rank/(nprocx*nprocy); - jproc = (rank-nprocx*nprocy*kproc)/nprocx; - iproc = rank-nprocx*nprocy*kproc-nprocz*jproc; + int kproc = rank/(nprocx*nprocy); + int jproc = (rank-nprocx*nprocy*kproc)/nprocx; + int iproc = rank-nprocx*nprocy*kproc-nprocz*jproc; if (rank == 0) { printf("i,j,k proc=%d %d %d \n",iproc,jproc,kproc); @@ -125,9 +88,8 @@ int main(int argc, char **argv) } double iVol_global = 1.0/Nx/Ny/Nz/nprocx/nprocy/nprocz; - int BoundaryCondition=0; - Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); + Domain Dm(db); Nx += 2; diff --git a/tests/TestSegDist.cpp b/tests/TestSegDist.cpp index 11e8c2a7..be7d00a4 100644 --- a/tests/TestSegDist.cpp +++ b/tests/TestSegDist.cpp @@ -14,6 +14,20 @@ #include "analysis/eikonal.h" +std::shared_ptr loadInputs( int nprocs ) +{ + INSIST(nprocs==8, "TestSegDist: Number of MPI processes must be equal to 8"); + auto db = std::make_shared( ); + db->putScalar( "BC", 0 ); + db->putVector( "nproc", { 2, 2, 2 } ); + db->putVector( "n", { 50, 50, 50 } ); + db->putScalar( "nspheres", 0 ); + db->putVector( "L", { 1, 1, 1 } ); + return db; +} + + +//*************************************************************************************** int main(int argc, char **argv) { // Initialize MPI @@ -25,39 +39,35 @@ int main(int argc, char **argv) { int i,j,k,n,nn; int iproc,jproc,kproc; - int nx,ny,nz; - int Nx, Ny, Nz, N; - int nprocx, nprocy, nprocz, nspheres; - double Lx, Ly, Lz; - Nx = Ny = Nz = 50; - nx = ny = nz = 50; - N = Nx*Ny*Nz; - nprocx=nprocy=nprocz=2; - Lx = Ly = Lz = 1.0; - int BC=0; - if (nprocs != 8){ - ERROR("TestSegDist: Number of MPI processes must be equal to 8"); - } - if (nprocx !=2 || nprocz !=2 || nprocy !=2 ){ - ERROR("TestSegDist: MPI process grid must be 2x2x2"); - } + + // Load inputs + auto db = loadInputs( nprocs ); + int Nx = db->getVector( "n" )[0]; + int Ny = db->getVector( "n" )[1]; + int Nz = db->getVector( "n" )[2]; + int nprocx = db->getVector( "nproc" )[0]; + int nprocy = db->getVector( "nproc" )[1]; + int nprocz = db->getVector( "nproc" )[2]; + // Get the rank info - Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC); - for (k=0;k loadInputs( int nprocs ) +{ + auto db = std::make_shared( "Domain.in" ); + db->putScalar( "BC", 0 ); + db->putVector( "nproc", { 1, 1, 1 } ); + db->putVector( "n", { 100, 100, 100 } ); + db->putScalar( "nspheres", 1 ); + db->putVector( "L", { 1, 1, 1 } ); + return db; +} int main(int argc, char **argv) @@ -34,65 +41,33 @@ int main(int argc, char **argv) //....................................................................... // Reading the domain information file //....................................................................... - int nprocx, nprocy, nprocz, nx, ny, nz, nspheres; - double Lx, Ly, Lz; - int Nx,Ny,Nz; int i,j,k,n; - if (rank==0){ - /* ifstream domain("Domain.in"); - domain >> nprocx; - domain >> nprocy; - domain >> nprocz; - domain >> nx; - domain >> ny; - domain >> nz; - domain >> nspheres; - domain >> Lx; - domain >> Ly; - domain >> Lz; - */ - // Set the domain for single processor test - nprocx=nprocy=nprocz=1; - nx=ny=nz=100; - nspheres=1; - Lx=Ly=Lz=1; - } - MPI_Barrier(comm); - // Computational domain - MPI_Bcast(&nx,1,MPI_INT,0,comm); - MPI_Bcast(&ny,1,MPI_INT,0,comm); - MPI_Bcast(&nz,1,MPI_INT,0,comm); - MPI_Bcast(&nprocx,1,MPI_INT,0,comm); - MPI_Bcast(&nprocy,1,MPI_INT,0,comm); - MPI_Bcast(&nprocz,1,MPI_INT,0,comm); - MPI_Bcast(&nspheres,1,MPI_INT,0,comm); - MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm); - //................................................. - MPI_Barrier(comm); + // Load inputs + auto db = loadInputs( nprocs ); + int Nx = db->getVector( "n" )[0]; + int Ny = db->getVector( "n" )[1]; + int Nz = db->getVector( "n" )[2]; + int nprocx = db->getVector( "nproc" )[0]; + int nprocy = db->getVector( "nproc" )[1]; + int nprocz = db->getVector( "nproc" )[2]; + + if (rank==0){ + printf("********************************************************\n"); + printf("Sub-domain size = %i x %i x %i\n",Nx,Ny,Nz); + printf("********************************************************\n"); + } + - // Check that the number of processors >= the number of ranks - if ( rank==0 ) { - printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz); - printf("Number of MPI ranks used: %i \n", nprocs); - printf("Full domain size: %i x %i x %i \n",nx*nprocx,ny*nprocy,nz*nprocz); - } - if ( nprocs < nprocx*nprocy*nprocz ) - ERROR("Insufficient number of processors"); - // Set up the domain - int BC=0; // Get the rank info - Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC); + Domain Dm(db); // const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); TwoPhase Averages(Dm); - int N = (nx+2)*(ny+2)*(nz+2); - Nx = nx+2; - Ny = ny+2; - Nz = nz+2; - if (rank == 0) cout << "Domain set." << endl; + Nx += 2; + Ny += 2; + Nz += 2; + int N = Nx*Ny*Nz; //....................................................................... for ( k=1;k(1e-7*fabs(A+B),1e-20); } @@ -85,7 +85,7 @@ void testWriter( const std::string& format, std::vector& mes // Get a list of the timesteps PROFILE_START(format+"-read-timesteps"); - std::vector timesteps = IO::readTimesteps( path + "/" + summary_name ); + auto timesteps = IO::readTimesteps( path + "/" + summary_name ); PROFILE_STOP(format+"-read-timesteps"); if ( timesteps.size()==2 ) ut.passes(format+": Corrent number of timesteps"); @@ -117,7 +117,7 @@ void testWriter( const std::string& format, std::vector& mes pass = true; for (size_t k=0; k mesh = IO::getMesh(path,timestep,database,k); + auto mesh = IO::getMesh(path,timestep,database,k); PROFILE_STOP(format+"-read-getMesh"); if ( mesh.get()==NULL ) { printf("Failed to load %s\n",database.name.c_str()); @@ -126,7 +126,7 @@ void testWriter( const std::string& format, std::vector& mes } if ( database.name=="pointmesh" ) { // Check the pointmesh - std::shared_ptr pmesh = IO::getPointList(mesh); + auto pmesh = IO::getPointList(mesh); if ( pmesh.get()==NULL ) { pass = false; break; @@ -138,8 +138,8 @@ void testWriter( const std::string& format, std::vector& mes } if ( database.name=="trimesh" || database.name=="trilist" ) { // Check the trimesh/trilist - std::shared_ptr mesh1 = IO::getTriMesh(mesh); - std::shared_ptr mesh2 = IO::getTriList(mesh); + auto mesh1 = IO::getTriMesh(mesh); + auto mesh2 = IO::getTriList(mesh); if ( mesh1.get()==NULL || mesh2.get()==NULL ) { pass = false; break; @@ -198,8 +198,7 @@ void testWriter( const std::string& format, std::vector& mes auto mesh = IO::getMesh(path,timestep,database,k); for (size_t v=0; vvars.size(); v++) { PROFILE_START(format+"-read-getVariable"); - std::shared_ptr variable = - IO::getVariable(path,timestep,database,k,mesh0->vars[v]->name); + auto variable = IO::getVariable(path,timestep,database,k,mesh0->vars[v]->name); if ( format=="new" ) IO::reformatVariable( *mesh, *variable ); PROFILE_STOP(format+"-read-getVariable"); @@ -251,19 +250,19 @@ int main(int argc, char **argv) }; // Create the meshes - std::shared_ptr set1( new IO::PointList(N_points) ); + auto set1 = std::make_shared(N_points); for (int i=0; ipoints[i].x = x[i]; set1->points[i].y = y[i]; set1->points[i].z = z[i]; } - std::shared_ptr trimesh( new IO::TriMesh(N_tri,set1) ); + auto trimesh = std::make_shared(N_tri,set1); for (int i=0; iA[i] = tri[i][0]; trimesh->B[i] = tri[i][1]; trimesh->C[i] = tri[i][2]; } - std::shared_ptr trilist( new IO::TriList(*trimesh) ); + auto trilist = std::make_shared(*trimesh); for (int i=0; i domain( new IO::DomainMesh(rank_data,6,7,8,1.0,1.0,1.0) ); + auto domain = std::make_shared(rank_data,6,7,8,1.0,1.0,1.0); // Create the variables - const auto NodeVar = IO::VariableType::NodeVariable; - const auto VolVar = IO::VariableType::VolumeVariable; - std::shared_ptr set_node_mag( new IO::Variable(1,NodeVar,"Node_set_mag") ); - std::shared_ptr set_node_vec( new IO::Variable(3,NodeVar,"Node_set_vec") ); - std::shared_ptr list_node_mag( new IO::Variable(1,NodeVar,"Node_list_mag") ); - std::shared_ptr list_node_vec( new IO::Variable(3,NodeVar,"Node_list_vec") ); - std::shared_ptr point_node_mag( new IO::Variable(1,NodeVar,"Node_point_mag") ); - std::shared_ptr point_node_vec( new IO::Variable(3,NodeVar,"Node_point_vec") ); - std::shared_ptr domain_node_mag( new IO::Variable(1,NodeVar,"Node_domain_mag") ); - std::shared_ptr domain_node_vec( new IO::Variable(3,NodeVar,"Node_domain_vec") ); - std::shared_ptr set_cell_mag( new IO::Variable(1,VolVar,"Cell_set_mag") ); - std::shared_ptr set_cell_vec( new IO::Variable(3,VolVar,"Cell_set_vec") ); - std::shared_ptr list_cell_mag( new IO::Variable(1,VolVar,"Cell_list_mag") ); - std::shared_ptr list_cell_vec( new IO::Variable(3,VolVar,"Cell_list_vec") ); - std::shared_ptr domain_cell_mag( new IO::Variable(1,VolVar,"Cell_domain_mag") ); - std::shared_ptr domain_cell_vec( new IO::Variable(3,VolVar,"Cell_domain_vec") ); + const auto NodeVar = IO::VariableType::NodeVariable; + const auto VolVar = IO::VariableType::VolumeVariable; + auto set_node_mag = std::make_shared(1,NodeVar,"Node_set_mag"); + auto set_node_vec = std::make_shared(3,NodeVar,"Node_set_vec"); + auto list_node_mag = std::make_shared(1,NodeVar,"Node_list_mag"); + auto list_node_vec = std::make_shared(3,NodeVar,"Node_list_vec"); + auto point_node_mag = std::make_shared(1,NodeVar,"Node_point_mag"); + auto point_node_vec = std::make_shared(3,NodeVar,"Node_point_vec"); + auto domain_node_mag = std::make_shared(1,NodeVar,"Node_domain_mag"); + auto domain_node_vec = std::make_shared(3,NodeVar,"Node_domain_vec"); + auto set_cell_mag = std::make_shared(1,VolVar,"Cell_set_mag"); + auto set_cell_vec = std::make_shared(3,VolVar,"Cell_set_vec"); + auto list_cell_mag = std::make_shared(1,VolVar,"Cell_list_mag"); + auto list_cell_vec = std::make_shared(3,VolVar,"Cell_list_vec"); + auto domain_cell_mag = std::make_shared(1,VolVar,"Cell_domain_mag"); + auto domain_cell_vec = std::make_shared(3,VolVar,"Cell_domain_vec"); point_node_mag->data.resize( N_points ); point_node_vec->data.resize( N_points, 3 ); for (int i=0; i loadInputs( ) +{ + auto db = std::make_shared( "Domain.in" ); + const int dim = 50; + db->putScalar( "BC", 0 ); + return db; +} + + +//*************************************************************************************** int main(int argc, char **argv) { //***************************************** @@ -24,8 +35,6 @@ int main(int argc, char **argv) MPI_Comm_size(comm,&nprocs); { // parallel domain size (# of sub-domains) - int nprocx,nprocy,nprocz; - int iproc,jproc,kproc; int sendtag,recvtag; //***************************************** // MPI ranks for all 18 neighbors @@ -54,55 +63,24 @@ int main(int argc, char **argv) // Variables that specify the computational domain string FILENAME; - int Nx,Ny,Nz; // local sub-domain size - int nspheres; // number of spheres in the packing - double Lx,Ly,Lz; // Domain length int i,j,k,n; // pmmc threshold values - if (rank==0){ - //....................................................................... - // Reading the domain information file - //....................................................................... - ifstream domain("Domain.in"); - domain >> nprocx; - domain >> nprocy; - domain >> nprocz; - domain >> Nx; - domain >> Ny; - domain >> Nz; - domain >> nspheres; - domain >> Lx; - domain >> Ly; - domain >> Lz; - //....................................................................... - } - // ************************************************************** - // Broadcast simulation parameters from rank 0 to all other procs - MPI_Barrier(comm); - // Computational domain - MPI_Bcast(&Nx,1,MPI_INT,0,comm); - MPI_Bcast(&Ny,1,MPI_INT,0,comm); - MPI_Bcast(&Nz,1,MPI_INT,0,comm); - MPI_Bcast(&nprocx,1,MPI_INT,0,comm); - MPI_Bcast(&nprocy,1,MPI_INT,0,comm); - MPI_Bcast(&nprocz,1,MPI_INT,0,comm); - MPI_Bcast(&nspheres,1,MPI_INT,0,comm); - MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm); - //................................................. - MPI_Barrier(comm); // ************************************************************** - if (nprocs != nprocx*nprocy*nprocz){ - printf("nprocx = %i \n",nprocx); - printf("nprocy = %i \n",nprocy); - printf("nprocz = %i \n",nprocz); - INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!"); - } - + // Load inputs + auto db = loadInputs( ); + int Nx = db->getVector( "n" )[0]; + int Ny = db->getVector( "n" )[1]; + int Nz = db->getVector( "n" )[2]; + int nprocx = db->getVector( "nproc" )[0]; + int nprocy = db->getVector( "nproc" )[1]; + int nprocz = db->getVector( "nproc" )[2]; + int kproc = rank/(nprocx*nprocy); + int jproc = (rank-nprocx*nprocy*kproc)/nprocx; + int iproc = rank-nprocx*nprocy*kproc-nprocx*jproc; + if (rank==0){ printf("********************************************************\n"); printf("Sub-domain size = %i x %i x %i\n",Nz,Nz,Nz); @@ -111,7 +89,7 @@ int main(int argc, char **argv) } // Initialized domain and averaging framework for Two-Phase Flow - Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC); + Domain Dm(db); Dm.CommInit(comm); TwoPhase Averages(Dm); diff --git a/tests/lbpm_dfh_simulator.cpp b/tests/lbpm_dfh_simulator.cpp index 1ad4c245..d5568b13 100644 --- a/tests/lbpm_dfh_simulator.cpp +++ b/tests/lbpm_dfh_simulator.cpp @@ -96,7 +96,7 @@ int main(int argc, char **argv) int RESTART_INTERVAL=20000; //int ANALYSIS_)INTERVAL=1000; int BLOB_ANALYSIS_INTERVAL=1000; - int timestep = 0; + int timestep = 6; if (rank==0){ //............................................................. @@ -492,7 +492,7 @@ int main(int argc, char **argv) ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int)*Np); ScaLBL_DeviceBarrier(); delete [] TmpMap; - + // Compute the solid interaction potential and copy result to device if (rank==0) printf("Computing solid interaction potential \n"); double *Tmp; @@ -523,14 +523,14 @@ int main(int argc, char **argv) for (int kk=0; kk<5; kk++){ for (int jj=0; jj<5; jj++){ for (int ii=0; ii<5; ii++){ - + int index = kk*25+jj*5+ii; double distval= Dst[index]; int idi=i+ii-2; int idj=j+jj-2; int idk=k+kk-2; - + if (idi < 0) idi=0; if (idj < 0) idj=0; if (idk < 0) idk=0; @@ -540,16 +540,16 @@ int main(int argc, char **argv) int nn = idk*Nx*Ny + idj*Nx + idi; if (!(Mask.id[nn] > 0)){ - double vec_x = double(ii-2); - double vec_y = double(jj-2); - double vec_z = double(kk-2); - - double ALPHA=PhaseLabel[nn]; - double GAMMA=-2.f; - if (distval > 2.f) ALPHA=0.f; // symmetric cutoff distance - phi_x += ALPHA*exp(GAMMA*distval)*vec_x/distval; - phi_y += ALPHA*exp(GAMMA*distval)*vec_y/distval; - phi_z += ALPHA*exp(GAMMA*distval)*vec_z/distval; + double vec_x = double(ii-2); + double vec_y = double(jj-2); + double vec_z = double(kk-2); + + double ALPHA=PhaseLabel[nn]; + double GAMMA=-2.f; + if (distval > 2.f) ALPHA=0.f; // symmetric cutoff distance + phi_x += ALPHA*exp(GAMMA*distval)*vec_x/distval; + phi_y += ALPHA*exp(GAMMA*distval)*vec_y/distval; + phi_z += ALPHA*exp(GAMMA*distval)*vec_z/distval; } } } @@ -557,17 +557,17 @@ int main(int argc, char **argv) Tmp[idx] = phi_x; Tmp[idx+Np] = phi_y; Tmp[idx+2*Np] = phi_z; - + /* double d = Averages->SDs(n); double dx = Averages->SDs_x(n); double dy = Averages->SDs_y(n); double dz = Averages->SDs_z(n); double value=cns*exp(-bns*fabs(d))-cws*exp(-bns*fabs(d)); - + Tmp[idx] = value*dx; Tmp[idx+Np] = value*dy; Tmp[idx+2*Np] = value*dz; - */ + */ } } } @@ -642,7 +642,7 @@ int main(int argc, char **argv) // Copy the phase from the GPU -> CPU //........................................................................... ScaLBL_DeviceBarrier(); - ScaLBL_Comm.RegularLayout(Map,Phi,Averages->Phase); + ScaLBL_CopyToHost(Averages->Phase.data(),Phi,Np*sizeof(double)); ScaLBL_Comm.RegularLayout(Map,Pressure,Averages->Press); ScaLBL_Comm.RegularLayout(Map,&Velocity[0],Averages->Vel_x); ScaLBL_Comm.RegularLayout(Map,&Velocity[Np],Averages->Vel_y); @@ -751,7 +751,7 @@ int main(int argc, char **argv) PROFILE_STOP("Update"); // Run the analysis - analysis.run( timestep, *Averages, Phi, Pressure, Velocity, fq, Den ); + // analysis.run( timestep, *Averages, Phi, Pressure, Velocity, fq, Den ); } analysis.finish(); @@ -782,8 +782,7 @@ int main(int argc, char **argv) FILE *OUTFILE; sprintf(LocalRankFilename,"Phase.%05i.raw",rank); OUTFILE = fopen(LocalRankFilename,"wb"); - //fwrite(PhaseField.data(),8,N,OUTFILE); - fwrite(Averages->Phase.data(),8,N,OUTFILE); + fwrite(PhaseField.data(),8,N,OUTFILE); fclose(OUTFILE); PROFILE_STOP("Main"); diff --git a/tests/lbpm_segmented_decomp.cpp b/tests/lbpm_segmented_decomp.cpp index c6804988..05770f25 100644 --- a/tests/lbpm_segmented_decomp.cpp +++ b/tests/lbpm_segmented_decomp.cpp @@ -274,7 +274,7 @@ int main(int argc, char **argv) n=k*nx*ny+j*nx+i; total++; Dm.id[n] = newIDs[n]; - if (!(Dm.id[n] > 0)){ + if (Dm.id[n] == 0){ count++; } } @@ -307,7 +307,7 @@ int main(int argc, char **argv) } } float porosity = float(totalGlobal-countGlobal)/totalGlobal; - //printf("Original Porosity=%f\n",porosity); + printf("Original Porosity=%f\n",porosity); } count = 0; diff --git a/tests/lbpm_segmented_pp.cpp b/tests/lbpm_segmented_pp.cpp index 605fe331..3bdbf094 100644 --- a/tests/lbpm_segmented_pp.cpp +++ b/tests/lbpm_segmented_pp.cpp @@ -228,40 +228,18 @@ int main(int argc, char **argv) // DoubleArray Distance(nx,ny,nz); // DoubleArray Phase(nx,ny,nz); + int count = 0; // Solve for the position of the solid phase for (k=0;k 0)) id[n] = 0; + if (Dm.id[n] == 0) id[n] = 0; else id[n] = 1; } } } - - - int count = 0; - int total = 0; - int totalGlobal = 0; - int countGlobal = 0; - for (k=1;k 0)){ - count++; - } - } - } - } - MPI_Allreduce(&total,&totalGlobal,1,MPI_INT,MPI_SUM,comm); - MPI_Allreduce(&count,&countGlobal,1,MPI_INT,MPI_SUM,comm); - float poros = float(totalGlobal-countGlobal)/totalGlobal; - if (rank==0) printf("Porosity=%f\n",poros); - - // Initialize the signed distance function for (k=0;k