OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Adding Database

Browse Source
This commit is contained in:
Mark Berrill 2018-05-15 10:01:14 -04:00
parent cd67724ab0
commit d3e91bc829
25 changed files with 2412 additions and 1275 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
CMakeLists.txt
View File

@ -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 )

13
IO/Writer.cpp
View File

@ -121,7 +121,7 @@ static std::vector<IO::MeshDatabase> writeMeshesOrigFormat( const std::vector<IO
static IO::MeshDatabase getDatabase( const std::string& filename, const IO::MeshDataStruct& mesh, int format )
{
int rank = MPI_WORLD_RANK();
char domainname[10];
char domainname[100];
sprintf(domainname,"%s_%05i",mesh.meshName.c_str(),rank);
// Create the MeshDatabase
IO::MeshDatabase database;
@ -310,12 +310,11 @@ static void writeSiloDomainMesh( DBfile *fid, const IO::MeshDataStruct& meshData
info.jy*mesh.Ly/info.ny, (info.jy+1)*mesh.Ly/info.ny,
info.kz*mesh.Lz/info.nz, (info.kz+1)*mesh.Lz/info.nz };
std::array<int,3> N = { mesh.nx, mesh.ny, mesh.nz };
std::array<int,3> 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<int>( fid, meshname+"_rankinfo", { mesh.rank, mesh.nprocx, mesh.nprocy, mesh.nprocz } );
for (size_t i=0; i<meshData.vars.size(); i++) {
const IO::Variable& var = *meshData.vars[i];
const auto& var = *meshData.vars[i];
auto type = static_cast<silo::VariableType>( 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<IO::MeshDatabase> writeMeshesSilo(
sprintf(fullpath,"%s/%s",path.c_str(),filename);
auto fid = silo::open( fullpath, silo::CREATE );
for (size_t i=0; i<meshData.size(); i++) {
std::shared_ptr<IO::Mesh> mesh = meshData[i].mesh;
auto mesh = meshData[i].mesh;
meshes_written.push_back( write_domain_silo(fid,filename,meshData[i],format) );
}
silo::close( fid );

29
analysis/eikonal.hpp → analysis/eikonal.cpp
View File

@ -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<float> &Distance, const Array<char> &ID, const Domain &Dm, const int timesteps)
float Eikonal3D( Array<float> &Distance, const Array<char> &ID, const Domain &Dm, const int timesteps)
{
PROFILE_START("Eikonal3D");
@ -327,7 +323,7 @@ inline float Eikonal3D( Array<float> &Distance, const Array<char> &ID, const Dom
/******************************************************************
* A fast distance calculation *
******************************************************************/
inline bool CalcDist3DIteration( Array<float> &Distance, const Domain &Dm )
bool CalcDist3DIteration( Array<float> &Distance, const Domain &Dm )
{
const float sq2 = sqrt(2.0f);
const float sq3 = sqrt(3.0f);
@ -367,7 +363,7 @@ inline bool CalcDist3DIteration( Array<float> &Distance, const Domain &Dm )
}
return changed;
}
inline void CalcDist3D( Array<float> &Distance, const Array<char> &ID, const Domain &Dm )
void CalcDist3D( Array<float> &Distance, const Array<char> &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<float> &Distance, const Array<char> &ID, const Dom
/******************************************************************
* A fast distance calculation *
******************************************************************/
inline void CalcDistMultiLevelHelper( Array<float> &Distance, const Domain &Dm )
void CalcDistMultiLevelHelper( Array<float> &Distance, const Domain &Dm )
{
size_t ratio = 4;
std::function<float(const Array<float>&)> coarsen = [ratio]( const Array<float>& data )
@ -440,7 +436,10 @@ inline void CalcDistMultiLevelHelper( Array<float> &Distance, const Domain &Dm )
// Coarsen
Array<float> 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<int>( "n" );
db->putVector<int>( "n", { n[0]/ratio, n[1]/ratio, n[2]/ratio } );
Domain Dm2(db);
Dm2.CommInit(Dm.Comm);
fillHalo<float> 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<float> &Distance, const Domain &Dm )
}
}
}
inline void CalcDistMultiLevel( Array<float> &Distance, const Array<char> &ID, const Domain &Dm )
void CalcDistMultiLevel( Array<float> &Distance, const Array<char> &ID, const Domain &Dm )
{
PROFILE_START("Calc Distance Multilevel");
int Nx = Dm.Nx-2;
@ -515,5 +514,3 @@ inline void CalcDistMultiLevel( Array<float> &Distance, const Array<char> &ID, c
Distance = imfilter::imfilter_separable<float>( Distance, {H,H,H}, BC );
PROFILE_STOP("Calc Distance Multilevel");
}
#endif

6
analysis/eikonal.h
View File

@ -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<float> &Distance, const Array<char> &ID, const Domain &Dm, const int timesteps);
@ -52,7 +53,4 @@ void CalcDist3D( Array<float> &Distance, const Array<char> &ID, const Domain &Dm
void CalcDistMultiLevel( Array<float> &Distance, const Array<char> &ID, const Domain &Dm );
#include "analysis/eikonal.hpp"
#endif

164
cmake/macros.cmake
View File

@ -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} $<TARGET_OBJECTS:${tmp}> )
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#

597
common/Database.cpp Normal file
View File

@ -0,0 +1,597 @@
#include "common/Database.h"
#include "common/Utilities.h"
#include <algorithm>
#include <cstring>
#include <iomanip>
#include <sstream>
#include <string>
#include <tuple>
/********************************************************************
* 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<KeyData> Database::clone() const { return cloneDatabase(); }
std::shared_ptr<Database> Database::cloneDatabase() const
{
auto db = std::make_shared<Database>();
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<KeyData> 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<const KeyData> Database::getData( const std::string& key ) const
{
return const_cast<Database*>( this )->getData( key );
}
bool Database::isDatabase( const std::string& key ) const
{
auto ptr = getData( key );
auto ptr2 = std::dynamic_pointer_cast<const Database>( ptr );
return ptr2 != nullptr;
}
std::shared_ptr<Database> Database::getDatabase( const std::string& key )
{
std::shared_ptr<KeyData> ptr = getData( key );
std::shared_ptr<Database> ptr2 = std::dynamic_pointer_cast<Database>( 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<const Database> Database::getDatabase( const std::string& key ) const
{
return const_cast<Database*>( this )->getDatabase( key );
}
std::vector<std::string> Database::getAllKeys() const
{
std::vector<std::string> 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<Database> db )
{
d_data[key] = std::move( db );
}
void Database::putData( const std::string& key, std::shared_ptr<KeyData> data )
{
d_data[key] = std::move( data );
}
/********************************************************************
* Is the data of the given type *
********************************************************************/
template<>
bool Database::isType<double>( const std::string& key ) const
{
auto type = getData( key )->type();
return type == "double";
}
template<>
bool Database::isType<float>( const std::string& key ) const
{
auto type = getData( key )->type();
return type == "double";
}
template<>
bool Database::isType<int>( const std::string& key ) const
{
bool pass = true;
auto type = getData( key )->type();
if ( type == "double" ) {
auto data = getVector<double>( key );
for ( auto tmp : data )
pass = pass && static_cast<double>( static_cast<int>( tmp ) ) == tmp;
} else {
pass = false;
}
return pass;
}
template<>
bool Database::isType<std::string>( const std::string& key ) const
{
auto type = getData( key )->type();
return type == "string";
}
template<>
bool Database::isType<bool>( const std::string& key ) const
{
auto type = getData( key )->type();
return type == "bool";
}
/********************************************************************
* Get a vector *
********************************************************************/
template<>
std::vector<std::string> Database::getVector<std::string>(
const std::string& key, const Units& ) const
{
std::shared_ptr<const KeyData> ptr = getData( key );
if ( std::dynamic_pointer_cast<const EmptyKeyData>( ptr ) )
return std::vector<std::string>();
const auto* ptr2 = dynamic_cast<const KeyDataString*>( ptr.get() );
if ( ptr2 == nullptr ) {
ERROR( "Key '" + key + "' is not a string" );
}
return ptr2->d_data;
}
template<>
std::vector<bool> Database::getVector<bool>( const std::string& key, const Units& ) const
{
std::shared_ptr<const KeyData> ptr = getData( key );
if ( std::dynamic_pointer_cast<const EmptyKeyData>( ptr ) )
return std::vector<bool>();
const auto* ptr2 = dynamic_cast<const KeyDataBool*>( ptr.get() );
if ( ptr2 == nullptr ) {
ERROR( "Key '" + key + "' is not a bool" );
}
return ptr2->d_data;
}
template<class TYPE>
std::vector<TYPE> Database::getVector( const std::string& key, const Units& unit ) const
{
std::shared_ptr<const KeyData> ptr = getData( key );
if ( std::dynamic_pointer_cast<const EmptyKeyData>( ptr ) )
return std::vector<TYPE>();
std::vector<TYPE> data;
if ( std::dynamic_pointer_cast<const KeyDataDouble>( ptr ) ) {
const auto* ptr2 = dynamic_cast<const KeyDataDouble*>( ptr.get() );
const std::vector<double>& 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<TYPE>( factor * data2[i] );
} else if ( std::dynamic_pointer_cast<const KeyDataString>( ptr ) ) {
ERROR( "Converting std::string to another type" );
} else if ( std::dynamic_pointer_cast<const KeyDataBool>( ptr ) ) {
ERROR( "Converting std::bool to another type" );
} else {
ERROR( "Unable to convert data format" );
}
return data;
}
/********************************************************************
* Put a vector *
********************************************************************/
template<>
void Database::putVector<std::string>(
const std::string& key, const std::vector<std::string>& data, const Units& )
{
std::shared_ptr<KeyDataString> ptr( new KeyDataString() );
ptr->d_data = data;
d_data[key] = ptr;
}
template<>
void Database::putVector<bool>(
const std::string& key, const std::vector<bool>& data, const Units& )
{
std::shared_ptr<KeyDataBool> ptr( new KeyDataBool() );
ptr->d_data = data;
d_data[key] = ptr;
}
template<class TYPE>
void Database::putVector( const std::string& key, const std::vector<TYPE>& data, const Units& unit )
{
std::shared_ptr<KeyDataDouble> 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<double>( 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<const Database*>( it.second.get() ) ) {
const auto* db = dynamic_cast<const Database*>( 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> Database::createFromString( const std::string& data )
{
std::shared_ptr<Database> 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<size_t, token_type> find_next_token( const char* buffer )
{
size_t i = 0;
while ( true ) {
if ( buffer[i] == '\n' || buffer[i] == '\r' ) {
return std::pair<size_t, token_type>( i + 1, token_type::newline );
} else if ( buffer[i] == 0 ) {
return std::pair<size_t, token_type>( i + 1, token_type::end );
} else if ( buffer[i] == '"' ) {
return std::pair<size_t, token_type>( i + 1, token_type::quote );
} else if ( buffer[i] == '=' ) {
return std::pair<size_t, token_type>( i + 1, token_type::equal );
} else if ( buffer[i] == '{' ) {
return std::pair<size_t, token_type>( i + 1, token_type::bracket );
} else if ( buffer[i] == '}' ) {
return std::pair<size_t, token_type>( i + 1, token_type::end_bracket );
} else if ( buffer[i] == '/' ) {
if ( buffer[i + 1] == '/' ) {
return std::pair<size_t, token_type>( i + 2, token_type::line_comment );
} else if ( buffer[i + 1] == '*' ) {
return std::pair<size_t, token_type>( i + 2, token_type::block_start );
}
} else if ( buffer[i] == '*' ) {
if ( buffer[i + 1] == '/' )
return std::pair<size_t, token_type>( i + 2, token_type::block_stop );
}
i++;
}
return std::pair<size_t, token_type>( 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<size_t, std::shared_ptr<KeyData>> 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<std::string> 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<KeyData> 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<size_t, std::shared_ptr<KeyData>>( 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<int>( 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<KeyData> 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> 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<double>::infinity() ) {
os << "inf";
} else if ( d_data[i] == -std::numeric_limits<double>::infinity() ) {
os << "-inf";
} else {
os << std::setprecision( 12 ) << d_data[i];
}
}
if ( !d_unit.isNull() )
os << " " << d_unit.str();
os << std::endl;
}
std::tuple<double, Units> 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<double>::infinity();
} else if ( tmp == "-inf" || tmp == "-infinity" ) {
data = -std::numeric_limits<double>::infinity();
} else if ( tmp == "nan" ) {
data = std::numeric_limits<double>::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<size_t>( pos - tmp.c_str() ) == tmp.size() + 1 )
ERROR( "Error reading value" );
}
return data;
}
/********************************************************************
* Instantiations *
********************************************************************/
template std::vector<char> Database::getVector<char>( const std::string&, const Units& ) const;
template std::vector<int> Database::getVector<int>( const std::string&, const Units& ) const;
template std::vector<size_t> Database::getVector<size_t>( const std::string&, const Units& ) const;
template std::vector<float> Database::getVector<float>( const std::string&, const Units& ) const;
template std::vector<double> Database::getVector<double>( const std::string&, const Units& ) const;
template void Database::putVector<char>(
const std::string&, const std::vector<char>&, const Units& );
template void Database::putVector<int>( const std::string&, const std::vector<int>&, const Units& );
template void Database::putVector<size_t>(
const std::string&, const std::vector<size_t>&, const Units& );
template void Database::putVector<float>(
const std::string&, const std::vector<float>&, const Units& );
template void Database::putVector<double>(
const std::string&, const std::vector<double>&, const Units& );
template bool Database::isType<int>( const std::string& ) const;
template bool Database::isType<float>( const std::string& ) const;
template bool Database::isType<double>( const std::string& ) const;
template bool Database::isType<std::string>( const std::string& ) const;

288
common/Database.h Normal file
View File

@ -0,0 +1,288 @@
#ifndef included_Database
#define included_Database
#include <iostream>
#include <map>
#include <memory>
#include <string>
#include <vector>
#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<KeyData> 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<Database> 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<KeyData> clone() const override;
//! Copy the data
std::shared_ptr<Database> 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<std::string> 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<class TYPE>
inline TYPE getScalar( const std::string& key, const Units& unit = Units() ) const;
/// @copydoc Database::getScalar(const std::string&,const Units&) const
template<class TYPE>
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<class TYPE>
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<class TYPE>
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<class TYPE>
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<class TYPE>
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<class TYPE>
std::vector<TYPE> getVector( const std::string& key, const Units& unit = Units() ) const;
/// @copydoc Database::getVector(const std::string&,const Units&) const
template<class TYPE>
inline std::vector<TYPE> 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<class TYPE>
void putVector(
const std::string& key, const std::vector<TYPE>& data, const Units& unit = Units() );
/// @copydoc Database::putVector(const std::string&,const std::vector<TYPE>&,const Units&)
template<class TYPE>
inline void putVector(
const std::string& key, const std::vector<TYPE>& 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<KeyData> 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<const KeyData> 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<KeyData> 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<class TYPE>
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<Database> 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<const Database> 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<Database> 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<std::string, std::shared_ptr<KeyData>> d_data;
// Function to load a database from a buffer
static size_t loadDatabase( const char* buffer, Database& db );
};
#include "common/Database.hpp"
#endif

173
common/Database.hpp Normal file
View File

@ -0,0 +1,173 @@
#ifndef included_Database_hpp
#define included_Database_hpp
#include "common/Database.h"
#include "common/Utilities.h"
#include <tuple>
/********************************************************************
* Basic classes for primative data types *
********************************************************************/
class EmptyKeyData : public KeyData
{
public:
EmptyKeyData() {}
virtual ~EmptyKeyData() {}
virtual std::shared_ptr<KeyData> clone() const override
{
return std::make_shared<EmptyKeyData>();
}
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<double>& data, const Units& unit )
: d_data( data ), d_unit( unit )
{
}
virtual ~KeyDataDouble() {}
virtual std::shared_ptr<KeyData> clone() const override
{
return std::make_shared<KeyDataDouble>( 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<double, Units> read( const std::string& );
static double readValue( const std::string& );
public:
std::vector<double> d_data;
Units d_unit;
};
class KeyDataBool : public KeyData
{
public:
KeyDataBool() {}
explicit KeyDataBool( const std::vector<bool>& data ) : d_data( data ) {}
virtual ~KeyDataBool() {}
virtual std::shared_ptr<KeyData> clone() const override
{
return std::make_shared<KeyDataBool>( 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<bool> d_data;
};
class KeyDataString : public KeyData
{
public:
KeyDataString() {}
explicit KeyDataString( const std::vector<std::string>& data ) : d_data( data ) {}
virtual ~KeyDataString() {}
virtual std::shared_ptr<KeyData> clone() const override
{
return std::make_shared<KeyDataString>( 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<std::string> d_data;
};
/********************************************************************
* Get a vector *
********************************************************************/
template<class TYPE>
inline std::vector<TYPE> Database::getVector(
const std::string& key, const std::string& unit ) const
{
return getVector<TYPE>( key, Units( unit ) );
}
template<class TYPE>
inline void Database::putVector(
const std::string& key, const std::vector<TYPE>& data, const std::string& unit )
{
putVector<TYPE>( key, data, Units( unit ) );
}
/********************************************************************
* Get a scalar *
********************************************************************/
template<class TYPE>
inline TYPE Database::getScalar( const std::string& key, const Units& unit ) const
{
const std::vector<TYPE>& data = getVector<TYPE>( 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<class TYPE>
inline TYPE Database::getWithDefault(
const std::string& key, const TYPE& value, const Units& unit ) const
{
if ( !keyExists( key ) )
return value;
return getScalar<TYPE>( key, unit );
}
template<class TYPE>
inline void Database::putScalar( const std::string& key, const TYPE& data, const Units& unit )
{
putVector<TYPE>( key, std::vector<TYPE>( 1, data ), unit );
}
template<class TYPE>
inline TYPE Database::getScalar( const std::string& key, const std::string& unit ) const
{
return getScalar<TYPE>( key, Units( unit ) );
}
template<class TYPE>
inline TYPE Database::getWithDefault(
const std::string& key, const TYPE& value, const std::string& unit ) const
{
return getWithDefault<TYPE>( key, value, Units( unit ) );
}
template<class TYPE>
inline void Database::putScalar( const std::string& key, const TYPE& data, const std::string& unit )
{
putScalar<TYPE>( key, data, Units( unit ) );
}
template<class TYPE>
inline void putVector(
const std::string& key, const std::vector<TYPE>& data, const std::string& unit )
{
putVector<TYPE>( key, data, Units( unit ) );
}
#endif

631
common/Domain.cpp
View File

@ -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<Database>( );
db->putScalar<int>( "BC", BC );
db->putVector<int>( "nproc", { npx, npx, npx } );
db->putVector<int>( "n", { nx, ny, nz } );
db->putScalar<int>( "nspheres", 0 );
db->putVector<double>( "L", { lx, ly, lz } );
initialize( db );
}
Domain::Domain( std::shared_ptr<Database> 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<Database> db )
{
d_db = db;
auto nproc = d_db->getVector<int>("nproc");
auto n = d_db->getVector<int>("n");
auto L = d_db->getVector<double>("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<int>("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<Q; q++){
weights[q] = sqrt(1.0*(D3Q27[q][0]*D3Q27[q][0]) + 1.0*(D3Q27[q][1]*D3Q27[q][1]) + 1.0*(D3Q27[q][2]*D3Q27[q][2]));
}
int xdim,ydim,zdim;
xdim=Dm.Nx-2;
ydim=Dm.Ny-2;
zdim=Dm.Nz-2;
fillHalo<double> 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<Dm.Nz-1;k++){
for (j=1;j<Dm.Ny-1;j++){
for (i=1;i<Dm.Nx-1;i++){
int n = k*Dm.Nx*Dm.Ny + j*Dm.Nx + i;
//sign = Distance(i,j,k) / fabs(Distance(i,j,k));
sign = -1;
if (ID[n] == 1) sign = 1;
/*
if (!(i+1<Nx)) nx=0.5*Distance(i,j,k);
else nx=0.5*Distance(i+1,j,k);;
if (!(j+1<Ny)) ny=0.5*Distance(i,j,k);
else ny=0.5*Distance(i,j+1,k);
if (!(k+1<Nz)) nz=0.5*Distance(i,j,k);
else nz=0.5*Distance(i,j,k+1);
if (i<1) nx-=0.5*Distance(i,j,k);
else nx-=0.5*Distance(i-1,j,k);
if (j<1) ny-=0.5*Distance(i,j,k);
else ny-=0.5*Distance(i,j-1,k);
if (k<1) nz-=0.5*Distance(i,j,k);
else nz-=0.5*Distance(i,j,k-1);
*/
//............Compute the Gradient...................................
nx = 0.5*(Distance(i+1,j,k) - Distance(i-1,j,k));
ny = 0.5*(Distance(i,j+1,k) - Distance(i,j-1,k));
nz = 0.5*(Distance(i,j,k+1) - Distance(i,j,k-1));
W = 0.0; Dx = Dy = Dz = 0.0;
// also ignore places where the gradient is zero since this will not
// result in any local change to Distance
if (nx*nx+ny*ny+nz*nz > 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;n<N;n++){
ID[n]=1;
}
//............................................
//............................................
// .........Loop over the spheres.............
for (p=0;p<nspheres;p++){
// Get the sphere from the list, map to local min
cx = List_cx[p] - min_x;
cy = List_cy[p] - min_y;
cz = List_cz[p] - min_z;
r = List_rad[p];
// Check if
// Range for this sphere in global indexing
imin = int ((cx-r)/hx)-1;
imax = int ((cx+r)/hx)+1;
jmin = int ((cy-r)/hy)-1;
jmax = int ((cy+r)/hy)+1;
kmin = int ((cz-r)/hz)-1;
kmax = int ((cz+r)/hz)+1;
// Obviously we have to do something at the edges
if (imin<0) imin = 0;
if (imin>Nx) 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<imax;i++){
for (j=jmin;j<jmax;j++){
for (k=kmin;k<kmax;k++){
// Initialize ID value to 'fluid (=1)'
x = i*hx;
y = j*hy;
z = k*hz;
value = 1;
// if inside sphere, set to zero
if ( (cx-x)*(cx-x)+(cy-y)*(cy-y)+(cz-z)*(cz-z) < r*r){
value=0;
}
// get the position in the list
n = k*Nx*Ny+j*Nx+i;
if ( ID[n] != 0 ){
ID[n] = value;
}
}
}
}
}
}
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)
{
// Use sphere lists to determine which nodes are in porespace
// Write out binary file for nodes
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 distance;
//............................................
// 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-2)*nprocx-1);
hy = Ly/((Ny-2)*nprocy-1);
hz = Lz/((Nz-2)*nprocz-1);
//............................................
// Get maximum and minimum for this domain
// Halo is included !
min_x = double(iproc*(Nx-2)-1)*hx;
min_y = double(jproc*(Ny-2)-1)*hy;
min_z = double(kproc*(Nz-2)-1)*hz;
//............................................
//............................................
// Pre-initialize Distance
for (n=0;n<N;n++){
Distance[n]=100.0;
}
//............................................
//............................................
// .........Loop over the spheres.............
for (p=0;p<nspheres;p++){
// Get the sphere from the list, map to local min
cx = List_cx[p] - min_x;
cy = List_cy[p] - min_y;
cz = List_cz[p] - min_z;
r = List_rad[p];
// Check if
// Range for this sphere in global indexing
imin = int ((cx-2*r)/hx);
imax = int ((cx+2*r)/hx)+2;
jmin = int ((cy-2*r)/hy);
jmax = int ((cy+2*r)/hy)+2;
kmin = int ((cz-2*r)/hz);
kmax = int ((cz+2*r)/hz)+2;
// Obviously we have to do something at the edges
if (imin<0) imin = 0;
if (imin>Nx) 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<imax;i++){
for (j=jmin;j<jmax;j++){
for (k=kmin;k<kmax;k++){
// x,y,z is distance in physical units
x = i*hx;
y = j*hy;
z = k*hz;
// if inside sphere, set to zero
// get the position in the list
n = k*Nx*Ny+j*Nx+i;
// Compute the distance
distance = sqrt((cx-x)*(cx-x)+(cy-y)*(cy-y)+(cz-z)*(cz-z)) - r;
// Assign the minimum distance
if (distance < Distance[n]) Distance[n] = distance;
}
}
}
}
// Map the distance to lattice units
for (n=0; n<N; n++) Distance[n] = Distance[n]/hx;
}
void WriteLocalSolidID(char *FILENAME, char *ID, int N)
{
char value;
ofstream File(FILENAME,ios::binary);
for (int n=0; n<N; n++){
value = ID[n];
File.write((char*) &value, sizeof(value));
}
File.close();
}
void WriteLocalSolidDistance(char *FILENAME, double *Distance, int N)
{
double value;
ofstream File(FILENAME,ios::binary);
for (int n=0; n<N; n++){
value = Distance[n];
File.write((char*) &value, sizeof(value));
}
File.close();
}
void WriteCheckpoint(const char *FILENAME, const double *cDen, const double *cfq, int Np)
{
int q,n;
double value;
ofstream File(FILENAME,ios::binary);
for (n=0; n<Np; n++){
// Write the two density values
value = cDen[n];
File.write((char*) &value, sizeof(value));
value = cDen[Np+n];
File.write((char*) &value, sizeof(value));
// Write the even distributions
for (q=0; q<19; q++){
value = cfq[q*Np+n];
File.write((char*) &value, sizeof(value));
}
}
File.close();
}
void ReadCheckpoint(char *FILENAME, double *cDen, double *cfq, int Np)
{
int q=0, n=0;
double value=0;
ifstream File(FILENAME,ios::binary);
for (n=0; n<Np; n++){
// Write the two density values
File.read((char*) &value, sizeof(value));
cDen[n] = value;
File.read((char*) &value, sizeof(value));
cDen[Np+n] = value;
// Read the even distributions
for (q=0; q<19; q++){
File.read((char*) &value, sizeof(value));
cfq[q*Np+n] = value;
}
}
File.close();
}
void ReadBinaryFile(char *FILENAME, double *Data, int N)
{
int n;
double value;
ifstream File(FILENAME,ios::binary);
if (File.good()){
for (n=0; n<N; n++){
// Write the two density values
File.read((char*) &value, sizeof(value));
Data[n] = value;
}
}
else {
for (n=0; n<N; n++) Data[n] = 1.2e-34;
}
File.close();
}

496
common/Domain.h
View File

@ -3,39 +3,60 @@
// Created by James McClure
// Copyright 2008-2013
//#define _GLIBCXX_USE_CXX11_ABI 0
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <math.h>
#include <time.h>
#include <exception> // std::exception
#include <exception>
#include <stdexcept>
#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<Database> 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<Database> 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<const Database> getDatabase() const { return d_db; }
private:
void initialize( std::shared_ptr<Database> db );
std::shared_ptr<Database> 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<Q; q++){
weights[q] = sqrt(1.0*(D3Q27[q][0]*D3Q27[q][0]) + 1.0*(D3Q27[q][1]*D3Q27[q][1]) + 1.0*(D3Q27[q][2]*D3Q27[q][2]));
}
void ReadCheckpoint(char *FILENAME, double *cDen, double *cfq, int Np);
int xdim,ydim,zdim;
xdim=Dm.Nx-2;
ydim=Dm.Ny-2;
zdim=Dm.Nz-2;
fillHalo<double> 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<Dm.Nz-1;k++){
for (j=1;j<Dm.Ny-1;j++){
for (i=1;i<Dm.Nx-1;i++){
int n = k*Dm.Nx*Dm.Ny + j*Dm.Nx + i;
//sign = Distance(i,j,k) / fabs(Distance(i,j,k));
sign = -1;
if (ID[n] == 1) sign = 1;
/*
if (!(i+1<Nx)) nx=0.5*Distance(i,j,k);
else nx=0.5*Distance(i+1,j,k);;
if (!(j+1<Ny)) ny=0.5*Distance(i,j,k);
else ny=0.5*Distance(i,j+1,k);
if (!(k+1<Nz)) nz=0.5*Distance(i,j,k);
else nz=0.5*Distance(i,j,k+1);
if (i<1) nx-=0.5*Distance(i,j,k);
else nx-=0.5*Distance(i-1,j,k);
if (j<1) ny-=0.5*Distance(i,j,k);
else ny-=0.5*Distance(i,j-1,k);
if (k<1) nz-=0.5*Distance(i,j,k);
else nz-=0.5*Distance(i,j,k-1);
*/
//............Compute the Gradient...................................
nx = 0.5*(Distance(i+1,j,k) - Distance(i-1,j,k));
ny = 0.5*(Distance(i,j+1,k) - Distance(i,j-1,k));
nz = 0.5*(Distance(i,j,k+1) - Distance(i,j,k-1));
W = 0.0; Dx = Dy = Dz = 0.0;
// also ignore places where the gradient is zero since this will not
// result in any local change to Distance
if (nx*nx+ny*ny+nz*nz > 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;n<N;n++){
ID[n]=1;
}
//............................................
//............................................
// .........Loop over the spheres.............
for (p=0;p<nspheres;p++){
// Get the sphere from the list, map to local min
cx = List_cx[p] - min_x;
cy = List_cy[p] - min_y;
cz = List_cz[p] - min_z;
r = List_rad[p];
// Check if
// Range for this sphere in global indexing
imin = int ((cx-r)/hx)-1;
imax = int ((cx+r)/hx)+1;
jmin = int ((cy-r)/hy)-1;
jmax = int ((cy+r)/hy)+1;
kmin = int ((cz-r)/hz)-1;
kmax = int ((cz+r)/hz)+1;
// Obviously we have to do something at the edges
if (imin<0) imin = 0;
if (imin>Nx) 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<imax;i++){
for (j=jmin;j<jmax;j++){
for (k=kmin;k<kmax;k++){
// Initialize ID value to 'fluid (=1)'
x = i*hx;
y = j*hy;
z = k*hz;
value = 1;
// if inside sphere, set to zero
if ( (cx-x)*(cx-x)+(cy-y)*(cy-y)+(cz-z)*(cz-z) < r*r){
value=0;
}
// get the position in the list
n = k*Nx*Ny+j*Nx+i;
if ( ID[n] != 0 ){
ID[n] = value;
}
}
}
}
}
}
inline 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)
{
// Use sphere lists to determine which nodes are in porespace
// Write out binary file for nodes
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 distance;
//............................................
// 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-2)*nprocx-1);
hy = Ly/((Ny-2)*nprocy-1);
hz = Lz/((Nz-2)*nprocz-1);
//............................................
// Get maximum and minimum for this domain
// Halo is included !
min_x = double(iproc*(Nx-2)-1)*hx;
min_y = double(jproc*(Ny-2)-1)*hy;
min_z = double(kproc*(Nz-2)-1)*hz;
//............................................
//............................................
// Pre-initialize Distance
for (n=0;n<N;n++){
Distance[n]=100.0;
}
//............................................
//............................................
// .........Loop over the spheres.............
for (p=0;p<nspheres;p++){
// Get the sphere from the list, map to local min
cx = List_cx[p] - min_x;
cy = List_cy[p] - min_y;
cz = List_cz[p] - min_z;
r = List_rad[p];
// Check if
// Range for this sphere in global indexing
imin = int ((cx-2*r)/hx);
imax = int ((cx+2*r)/hx)+2;
jmin = int ((cy-2*r)/hy);
jmax = int ((cy+2*r)/hy)+2;
kmin = int ((cz-2*r)/hz);
kmax = int ((cz+2*r)/hz)+2;
// Obviously we have to do something at the edges
if (imin<0) imin = 0;
if (imin>Nx) 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<imax;i++){
for (j=jmin;j<jmax;j++){
for (k=kmin;k<kmax;k++){
// x,y,z is distance in physical units
x = i*hx;
y = j*hy;
z = k*hz;
// if inside sphere, set to zero
// get the position in the list
n = k*Nx*Ny+j*Nx+i;
// Compute the distance
distance = sqrt((cx-x)*(cx-x)+(cy-y)*(cy-y)+(cz-z)*(cz-z)) - r;
// Assign the minimum distance
if (distance < Distance[n]) Distance[n] = distance;
}
}
}
}
// Map the distance to lattice units
for (n=0; n<N; n++) Distance[n] = Distance[n]/hx;
}
inline void WriteLocalSolidID(char *FILENAME, char *ID, int N)
{
char value;
ofstream File(FILENAME,ios::binary);
for (int n=0; n<N; n++){
value = ID[n];
File.write((char*) &value, sizeof(value));
}
File.close();
}
inline void WriteLocalSolidDistance(char *FILENAME, double *Distance, int N)
{
double value;
ofstream File(FILENAME,ios::binary);
for (int n=0; n<N; n++){
value = Distance[n];
File.write((char*) &value, sizeof(value));
}
File.close();
}
inline void WriteCheckpoint(const char *FILENAME, const double *cDen, const double *cfq, int Np)
{
int q,n;
double value;
ofstream File(FILENAME,ios::binary);
for (n=0; n<Np; n++){
// Write the two density values
value = cDen[n];
File.write((char*) &value, sizeof(value));
value = cDen[Np+n];
File.write((char*) &value, sizeof(value));
// Write the even distributions
for (q=0; q<19; q++){
value = cfq[q*Np+n];
File.write((char*) &value, sizeof(value));
}
}
File.close();
}
inline void ReadCheckpoint(char *FILENAME, double *cDen, double *cfq, int Np)
{
int q=0, n=0;
double value=0;
ifstream File(FILENAME,ios::binary);
for (n=0; n<Np; n++){
// Write the two density values
File.read((char*) &value, sizeof(value));
cDen[n] = value;
File.read((char*) &value, sizeof(value));
cDen[Np+n] = value;
// Read the even distributions
for (q=0; q<19; q++){
File.read((char*) &value, sizeof(value));
cfq[q*Np+n] = value;
}
}
File.close();
}
inline void ReadBinaryFile(char *FILENAME, double *Data, int N)
{
int n;
double value;
ifstream File(FILENAME,ios::binary);
if (File.good()){
for (n=0; n<N; n++){
// Write the two density values
File.read((char*) &value, sizeof(value));
Data[n] = value;
}
}
else {
for (n=0; n<N; n++) Data[n] = 1.2e-34;
}
File.close();
}
void ReadBinaryFile(char *FILENAME, double *Data, int N);
#endif

226
common/Units.cpp Normal file
View File

@ -0,0 +1,226 @@
#include "common/Units.h"
#include "common/Utilities.h"
#include <algorithm>
#include <cmath>
#include <string>
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<char, 3> Units::str( UnitPrefix p ) noexcept
{
std::array<char, 3> str;
str[0] = d_prefixSymbol[static_cast<int8_t>( 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";
}

140
common/Units.h Normal file
View File

@ -0,0 +1,140 @@
#ifndef included_Units
#define included_Units
#include <array>
#include <iostream>
#include <map>
#include <memory>
#include <string>
#include <vector>
//! 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<int8_t>( x )];
}
//! Get a string representation of the units
std::string str() const;
//! Get a string representation for the prefix
static std::array<char, 3> 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

131
tests/TestColorSquareTube.cpp
View File

@ -9,7 +9,35 @@
#include "common/ScaLBL.h"
#include "common/MPI_Helpers.h"
using namespace std;
std::shared_ptr<Database> loadInputs( int nprocs )
{
auto db = std::make_shared<Database>( "Domain.in" );
const int dim = 50;
db->putScalar<int>( "BC", 0 );
if ( nprocs == 1 ){
db->putVector<int>( "nproc", { 1, 1, 1 } );
db->putVector<int>( "n", { 3, 1, 1 } );
db->putScalar<int>( "nspheres", 0 );
db->putVector<double>( "L", { 1, 1, 1 } );
} else if ( nprocs == 2 ) {
db->putVector<int>( "nproc", { 2, 1, 1 } );
db->putVector<int>( "n", { dim, dim, dim } );
db->putScalar<int>( "nspheres", 0 );
db->putVector<double>( "L", { 1, 1, 1 } );
} else if ( nprocs == 4 ) {
db->putVector<int>( "nproc", { 2, 2, 1 } );
db->putVector<int>( "n", { dim, dim, dim } );
db->putScalar<int>( "nspheres", 0 );
db->putVector<double>( "L", { 1, 1, 1 } );
} else if (nprocs==8){
db->putVector<int>( "nproc", { 2, 2, 2 } );
db->putVector<int>( "n", { dim, dim, dim } );
db->putScalar<int>( "nspheres", 0 );
db->putVector<double>( "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<int>( "n" )[0];
int Ny = db->getVector<int>( "n" )[1];
int Nz = db->getVector<int>( "n" )[2];
int nprocx = db->getVector<int>( "nproc" )[0];
int nprocy = db->getVector<int>( "nproc" )[1];
int nprocz = db->getVector<int>( "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<Nz-1;k++){

120
tests/TestCommD3Q19.cpp
View File

@ -11,6 +11,36 @@
using namespace std;
std::shared_ptr<Database> loadInputs( int nprocs )
{
auto db = std::make_shared<Database>( "Domain.in" );
const int dim = 50;
db->putScalar<int>( "BC", 0 );
if ( nprocs == 1 ){
db->putVector<int>( "nproc", { 1, 1, 1 } );
db->putVector<int>( "n", { 3, 1, 1 } );
db->putScalar<int>( "nspheres", 0 );
db->putVector<double>( "L", { 1, 1, 1 } );
} else if ( nprocs == 2 ) {
db->putVector<int>( "nproc", { 2, 1, 1 } );
db->putVector<int>( "n", { dim, dim, dim } );
db->putScalar<int>( "nspheres", 0 );
db->putVector<double>( "L", { 1, 1, 1 } );
} else if ( nprocs == 4 ) {
db->putVector<int>( "nproc", { 2, 2, 1 } );
db->putVector<int>( "n", { dim, dim, dim } );
db->putScalar<int>( "nspheres", 0 );
db->putVector<double>( "L", { 1, 1, 1 } );
} else if (nprocs==8){
db->putVector<int>( "nproc", { 2, 2, 2 } );
db->putVector<int>( "n", { dim, dim, dim } );
db->putScalar<int>( "nspheres", 0 );
db->putVector<double>( "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(&timestepMax,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<int>( "n" )[0];
int Ny = db->getVector<int>( "n" )[1];
int Nz = db->getVector<int>( "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,

82
tests/TestForceMoments.cpp
View File

@ -30,6 +30,22 @@ extern void PrintNeighborList(int * neighborList, int Np, int rank) {
}
}
std::shared_ptr<Database> loadInputs( int nprocs )
{
auto db = std::make_shared<Database>( "Domain.in" );
const int dim = 50;
db->putScalar<int>( "BC", 0 );
if ( nprocs == 1 ){
db->putVector<int>( "nproc", { 1, 1, 1 } );
db->putVector<int>( "n", { 3, 3, 3 } );
db->putScalar<int>( "nspheres", 0 );
db->putVector<double>( "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<int>( "n" )[0];
int Ny = db->getVector<int>( "n" )[1];
int Nz = db->getVector<int>( "n" )[2];
int nprocx = db->getVector<int>( "nproc" )[0];
int nprocy = db->getVector<int>( "nproc" )[1];
int nprocz = db->getVector<int>( "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;

31
tests/TestInterfaceSpeed.cpp
View File

@ -15,6 +15,20 @@
#define SPEED -1
#define PI 3.14159
std::shared_ptr<Database> loadInputs( int nprocs )
{
auto db = std::make_shared<Database>( );
const int dim = 50;
db->putScalar<int>( "BC", 0 );
db->putVector<int>( "nproc", { 1, 1, 1 } );
db->putVector<int>( "n", { N, N, N } );
db->putScalar<int>( "nspheres", 0 );
db->putVector<double>( "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<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = 1;
@ -42,9 +51,9 @@ int main (int argc, char *argv[])
TwoPhase Averages(Dm);
int timestep=0;
Nx = Dm.Nx;
Ny = Dm.Ny;
Nz = Dm.Nz;
int Nx = Dm.Nx;
int Ny = Dm.Ny;
int Nz = Dm.Nz;
double Cx,Cy,Cz;
double dist1,dist2;

118
tests/TestPoiseuille.cpp
View File

@ -9,7 +9,35 @@
#include "common/ScaLBL.h"
#include "common/MPI_Helpers.h"
using namespace std;
bool exists( const std::string& filename )
{
ifstream domain( filename );
return domain.good();
}
std::shared_ptr<Database> loadInputs( int nprocs )
{
const int dim = 12;
std::shared_ptr<Database> db;
if ( exists( "Domain.in" ) ) {
db = std::make_shared<Database>( "Domain.in" );
} else if (nprocs==1) {
db->putVector<int>( "nproc", { 1, 1, 1 } );
db->putVector<int>( "n", { dim, dim, dim } );
db->putScalar<int>( "nspheres", 0 );
db->putVector<double>( "L", { 1, 1, 1 } );
} else if (nprocs==2) {
db->putVector<int>( "nproc", { 1, 1, 1 } );
db->putVector<int>( "n", { dim/2, dim/2, dim/2 } );
db->putScalar<int>( "nspheres", 0 );
db->putVector<double>( "L", { 1, 1, 1 } );
}
db->putScalar<int>( "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<int>( "n" )[0];
int Ny = db->getVector<int>( "n" )[1];
int Nz = db->getVector<int>( "n" )[2];
int nprocx = db->getVector<int>( "nproc" )[0];
int nprocy = db->getVector<int>( "nproc" )[1];
int nprocz = db->getVector<int>( "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;

92
tests/TestPressVel.cpp
View File

@ -9,7 +9,21 @@
#include "common/ScaLBL.h"
#include "common/MPI_Helpers.h"
using namespace std;
std::shared_ptr<Database> loadInputs( int nprocs )
{
auto db = std::make_shared<Database>( "Domain.in" );
const int dim = 50;
db->putScalar<int>( "BC", 0 );
if ( nprocs == 1 ){
db->putVector<int>( "nproc", { 1, 1, 1 } );
db->putVector<int>( "n", { 4, 4, 4 } );
db->putScalar<int>( "nspheres", 0 );
db->putVector<double>( "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<int>( "n" )[0];
int Ny = db->getVector<int>( "n" )[1];
int Nz = db->getVector<int>( "n" )[2];
int nprocx = db->getVector<int>( "nproc" )[0];
int nprocy = db->getVector<int>( "nproc" )[1];
int nprocz = db->getVector<int>( "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;

56
tests/TestSegDist.cpp
View File

@ -14,6 +14,20 @@
#include "analysis/eikonal.h"
std::shared_ptr<Database> loadInputs( int nprocs )
{
INSIST(nprocs==8, "TestSegDist: Number of MPI processes must be equal to 8");
auto db = std::make_shared<Database>( );
db->putScalar<int>( "BC", 0 );
db->putVector<int>( "nproc", { 2, 2, 2 } );
db->putVector<int>( "n", { 50, 50, 50 } );
db->putScalar<int>( "nspheres", 0 );
db->putVector<double>( "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<int>( "n" )[0];
int Ny = db->getVector<int>( "n" )[1];
int Nz = db->getVector<int>( "n" )[2];
int nprocx = db->getVector<int>( "nproc" )[0];
int nprocy = db->getVector<int>( "nproc" )[1];
int nprocz = db->getVector<int>( "nproc" )[2];
// Get the rank info
Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
for (k=0;k<nz;k++){
for (j=0;j<ny;j++){
for (i=0;i<nx;i++){
n = k*nx*ny+j*nx+i;
Domain Dm(db);
for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
Dm.id[n] = 1;
}
}
}
Dm.CommInit(comm);
nx+=2; ny+=2; nz+=2;
N = nx*ny*nz;
int nx = Nx+2;
int ny = Ny+2;
int nz = Nz+2;
int N = nx*ny*nz;
int count = 0;
char *id;

85
tests/TestTorus.cpp
View File

@ -9,10 +9,17 @@
#include "analysis/analysis.h"
#include "analysis/TwoPhase.h"
//#include "Domain.h"
using namespace std;
std::shared_ptr<Database> loadInputs( int nprocs )
{
auto db = std::make_shared<Database>( "Domain.in" );
db->putScalar<int>( "BC", 0 );
db->putVector<int>( "nproc", { 1, 1, 1 } );
db->putVector<int>( "n", { 100, 100, 100 } );
db->putScalar<int>( "nspheres", 1 );
db->putVector<double>( "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<int>( "n" )[0];
int Ny = db->getVector<int>( "n" )[1];
int Nz = db->getVector<int>( "n" )[2];
int nprocx = db->getVector<int>( "nproc" )[0];
int nprocy = db->getVector<int>( "nproc" )[1];
int nprocz = db->getVector<int>( "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<Nz-1;k++){
for ( j=1;j<Ny-1;j++){

55
tests/TestWriter.cpp
View File

@ -21,7 +21,7 @@ inline bool approx_equal( const Point& A, const Point& B )
}
inline bool approx_equal( const double& A, const double& B )
{
return fabs(A-B)<=1e-7*fabs(A+B);
return fabs(A-B) <= std::max<double>(1e-7*fabs(A+B),1e-20);
}
@ -85,7 +85,7 @@ void testWriter( const std::string& format, std::vector<IO::MeshDataStruct>& mes
// Get a list of the timesteps
PROFILE_START(format+"-read-timesteps");
std::vector<std::string> 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<IO::MeshDataStruct>& mes
pass = true;
for (size_t k=0; k<database.domains.size(); k++) {
PROFILE_START(format+"-read-getMesh");
std::shared_ptr<IO::Mesh> 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<IO::MeshDataStruct>& mes
}
if ( database.name=="pointmesh" ) {
// Check the pointmesh
std::shared_ptr<IO::PointList> 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<IO::MeshDataStruct>& mes
}
if ( database.name=="trimesh" || database.name=="trilist" ) {
// Check the trimesh/trilist
std::shared_ptr<IO::TriMesh> mesh1 = IO::getTriMesh(mesh);
std::shared_ptr<IO::TriList> 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<IO::MeshDataStruct>& mes
auto mesh = IO::getMesh(path,timestep,database,k);
for (size_t v=0; v<mesh0->vars.size(); v++) {
PROFILE_START(format+"-read-getVariable");
std::shared_ptr<IO::Variable> 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<IO::PointList> set1( new IO::PointList(N_points) );
auto set1 = std::make_shared<IO::PointList>(N_points);
for (int i=0; i<N_points; i++) {
set1->points[i].x = x[i];
set1->points[i].y = y[i];
set1->points[i].z = z[i];
}
std::shared_ptr<IO::TriMesh> trimesh( new IO::TriMesh(N_tri,set1) );
auto trimesh = std::make_shared<IO::TriMesh>(N_tri,set1);
for (int i=0; i<N_tri; i++) {
trimesh->A[i] = tri[i][0];
trimesh->B[i] = tri[i][1];
trimesh->C[i] = tri[i][2];
}
std::shared_ptr<IO::TriList> trilist( new IO::TriList(*trimesh) );
auto trilist = std::make_shared<IO::TriList>(*trimesh);
for (int i=0; i<N_tri; i++) {
Point A(x[tri[i][0]],y[tri[i][0]],z[tri[i][0]]);
Point B(x[tri[i][1]],y[tri[i][1]],z[tri[i][1]]);
@ -275,25 +274,25 @@ int main(int argc, char **argv)
}
}
RankInfoStruct rank_data( rank, nprocs, 1, 1 );
std::shared_ptr<IO::DomainMesh> domain( new IO::DomainMesh(rank_data,6,7,8,1.0,1.0,1.0) );
auto domain = std::make_shared<IO::DomainMesh>(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<IO::Variable> set_node_mag( new IO::Variable(1,NodeVar,"Node_set_mag") );
std::shared_ptr<IO::Variable> set_node_vec( new IO::Variable(3,NodeVar,"Node_set_vec") );
std::shared_ptr<IO::Variable> list_node_mag( new IO::Variable(1,NodeVar,"Node_list_mag") );
std::shared_ptr<IO::Variable> list_node_vec( new IO::Variable(3,NodeVar,"Node_list_vec") );
std::shared_ptr<IO::Variable> point_node_mag( new IO::Variable(1,NodeVar,"Node_point_mag") );
std::shared_ptr<IO::Variable> point_node_vec( new IO::Variable(3,NodeVar,"Node_point_vec") );
std::shared_ptr<IO::Variable> domain_node_mag( new IO::Variable(1,NodeVar,"Node_domain_mag") );
std::shared_ptr<IO::Variable> domain_node_vec( new IO::Variable(3,NodeVar,"Node_domain_vec") );
std::shared_ptr<IO::Variable> set_cell_mag( new IO::Variable(1,VolVar,"Cell_set_mag") );
std::shared_ptr<IO::Variable> set_cell_vec( new IO::Variable(3,VolVar,"Cell_set_vec") );
std::shared_ptr<IO::Variable> list_cell_mag( new IO::Variable(1,VolVar,"Cell_list_mag") );
std::shared_ptr<IO::Variable> list_cell_vec( new IO::Variable(3,VolVar,"Cell_list_vec") );
std::shared_ptr<IO::Variable> domain_cell_mag( new IO::Variable(1,VolVar,"Cell_domain_mag") );
std::shared_ptr<IO::Variable> 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<IO::Variable>(1,NodeVar,"Node_set_mag");
auto set_node_vec = std::make_shared<IO::Variable>(3,NodeVar,"Node_set_vec");
auto list_node_mag = std::make_shared<IO::Variable>(1,NodeVar,"Node_list_mag");
auto list_node_vec = std::make_shared<IO::Variable>(3,NodeVar,"Node_list_vec");
auto point_node_mag = std::make_shared<IO::Variable>(1,NodeVar,"Node_point_mag");
auto point_node_vec = std::make_shared<IO::Variable>(3,NodeVar,"Node_point_vec");
auto domain_node_mag = std::make_shared<IO::Variable>(1,NodeVar,"Node_domain_mag");
auto domain_node_vec = std::make_shared<IO::Variable>(3,NodeVar,"Node_domain_vec");
auto set_cell_mag = std::make_shared<IO::Variable>(1,VolVar,"Cell_set_mag");
auto set_cell_vec = std::make_shared<IO::Variable>(3,VolVar,"Cell_set_vec");
auto list_cell_mag = std::make_shared<IO::Variable>(1,VolVar,"Cell_list_mag");
auto list_cell_vec = std::make_shared<IO::Variable>(3,VolVar,"Cell_list_vec");
auto domain_cell_mag = std::make_shared<IO::Variable>(1,VolVar,"Cell_domain_mag");
auto domain_cell_vec = std::make_shared<IO::Variable>(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<N_points; i++) {

70
tests/lbpm_captube_pp.cpp
View File

@ -11,6 +11,17 @@
#include "analysis/TwoPhase.h"
#include "common/MPI_Helpers.h"
std::shared_ptr<Database> loadInputs( )
{
auto db = std::make_shared<Database>( "Domain.in" );
const int dim = 50;
db->putScalar<int>( "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<int>( "n" )[0];
int Ny = db->getVector<int>( "n" )[1];
int Nz = db->getVector<int>( "n" )[2];
int nprocx = db->getVector<int>( "nproc" )[0];
int nprocy = db->getVector<int>( "nproc" )[1];
int nprocz = db->getVector<int>( "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);

41
tests/lbpm_dfh_simulator.cpp
View File

@ -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");

4
tests/lbpm_segmented_decomp.cpp
View File

@ -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;

26
tests/lbpm_segmented_pp.cpp
View File

@ -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<nz;k++){
for (j=0;j<ny;j++){
for (i=0;i<nx;i++){
n = k*nx*ny+j*nx+i;
// Initialize the solid phase
if (!(Dm.id[n] > 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<nz-1;k++){
for (j=1;j<ny-1;j++){
for (i=1;i<nx-1;i++){
n=k*nx*ny+j*nx+i;
total++;
if (!(id[n] > 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<nz;k++){
for (j=0;j<ny;j++){