Merge branch 'master' of github.com:JamesEMcClure/LBPM-WIA

This commit is contained in:
James E McClure
2016-02-28 17:39:41 -05:00
8 changed files with 559 additions and 12 deletions

View File

@@ -105,6 +105,7 @@ IF ( NOT ONLY_BUILD_DOCS )
CONFIGURE_MPI() # MPI must be before other libraries
CONFIGURE_CUDA()
CONFIGURE_MIC()
CONFIGURE_NETCDF()
CONFIGURE_LBPM()
CONFIGURE_TIMER( 0 "${${PROJ}_INSTALL_DIR}/null_timer" )
CONFIGURE_LINE_COVERAGE()

288
IO/netcdf.cpp Normal file
View File

@@ -0,0 +1,288 @@
#include "IO/netcdf.h"
#include "common/Utilities.h"
#include "ProfilerApp.h"
#ifdef USE_NETCDF
#include <netcdf.h>
#define CHECK_NC_ERR( ERR ) \
do { \
if ( ERR != NC_NOERR ) { \
std::string msg = "Error calling netcdf routine: "; \
msg += nc_strerror( ERR ); \
ERROR( msg ); \
} \
} while (0)
namespace netcdf {
// Function to reverse an array
template<class TYPE>
inline std::vector<TYPE> reverse( const std::vector<TYPE>& x )
{
std::vector<TYPE> y(x.size());
for (size_t i=0; i<x.size(); i++)
y[i] = x[x.size()-i-1];
return y;
}
/****************************************************
* Open/close a file *
****************************************************/
int open( const std::string& filename )
{
int fid = 0;
int err = nc_open( filename.c_str(), NC_NOWRITE, &fid );
CHECK_NC_ERR( err );
return fid;
}
void close( int fid )
{
int err = nc_close( fid );
if ( err != NC_NOERR )
ERROR("Error closing file");
}
/****************************************************
* Query basic properties *
****************************************************/
static std::vector<size_t> getDimVar( int fid, int varid )
{
int ndim = 0;
int err = nc_inq_varndims( fid, varid, &ndim );
CHECK_NC_ERR( err );
std::vector<size_t> dims(ndim,0);
int dimid[64] = {-1};
err = nc_inq_vardimid( fid, varid, dimid );
CHECK_NC_ERR( err );
for (int i=0; i<ndim; i++) {
err = nc_inq_dimlen( fid, dimid[i], &dims[i] );
CHECK_NC_ERR( err );
}
return dims;
}
static int getVarID( int fid, const std::string& var )
{
int id = -1;
int err = nc_inq_varid( fid, var.c_str(), &id );
CHECK_NC_ERR( err );
return id;
}
std::vector<size_t> getVarDim( int fid, const std::string& var )
{
return getDimVar( fid, getVarID( fid, var ) );
}
std::vector<size_t> getAttDim( int fid, const std::string& att )
{
std::vector<size_t> dim(1,0);
int err = nc_inq_attlen( fid, NC_GLOBAL, att.c_str(), dim.data() );
return dim;
}
std::vector<std::string> getVarNames( int fid )
{
int nvar;
int err = nc_inq( fid, NULL, &nvar, NULL, NULL );
CHECK_NC_ERR( err );
std::vector<std::string> vars(nvar);
for (int i=0; i<nvar; i++) {
char name[NC_MAX_NAME+1];
err = nc_inq_varname( fid, i, name );
CHECK_NC_ERR( err );
vars[i] = name;
}
return vars;
}
std::vector<std::string> getAttNames( int fid )
{
int natt;
int err = nc_inq( fid, NULL, NULL, &natt, NULL );
CHECK_NC_ERR( err );
std::vector<std::string> att(natt);
for (int i=0; i<natt; i++) {
char name[NC_MAX_NAME+1];
err = nc_inq_attname( fid, NC_GLOBAL, i, name );
CHECK_NC_ERR( err );
att[i] = name;
}
return att;
}
static inline VariableType convertType( nc_type type )
{
VariableType type2;
if ( type == NC_BYTE )
type2 = BYTE;
else if ( type == NC_CHAR )
type2 = STRING;
else if ( type == NC_SHORT )
type2 = SHORT;
else if ( type == NC_USHORT )
type2 = USHORT;
else if ( type == NC_INT )
type2 = INT;
else if ( type == NC_UINT )
type2 = UINT;
else if ( type == NC_INT64 )
type2 = INT64;
else if ( type == NC_UINT64 )
type2 = UINT64;
else if ( type == NC_FLOAT )
type2 = FLOAT;
else if ( type == NC_DOUBLE )
type2 = DOUBLE;
else
ERROR("Unknown type");
return type2;
}
VariableType getVarType( int fid, const std::string& var )
{
int varid = -1;
int err = nc_inq_varid( fid, var.c_str(), &varid );
CHECK_NC_ERR( err );
nc_type type;
err = nc_inq_vartype( fid, varid, &type );
CHECK_NC_ERR( err );
return convertType(type);
}
VariableType getAttType( int fid, const std::string& att )
{
nc_type type;
int err = nc_inq_atttype( fid, NC_GLOBAL, att.c_str(), &type );
CHECK_NC_ERR( err );
return convertType(type);
}
/****************************************************
* Read a variable *
****************************************************/
template<>
Array<unsigned short> getVar<unsigned short>( int fid, const std::string& var )
{
PROFILE_START("getVar<unsigned short>");
Array<unsigned short> x( reverse(getVarDim(fid,var)) );
int err = nc_get_var_ushort( fid, getVarID(fid,var), x.get() );
CHECK_NC_ERR( err );
PROFILE_STOP("getVar<unsigned short>");
return x;
}
template<>
Array<short> getVar<short>( int fid, const std::string& var )
{
PROFILE_START("getVar<short>");
Array<short> x( reverse(getVarDim(fid,var)) );
int err = nc_get_var_short( fid, getVarID(fid,var), x.get() );
CHECK_NC_ERR( err );
PROFILE_STOP("getVar<short>");
return x;
}
template<>
Array<unsigned int> getVar<unsigned int>( int fid, const std::string& var )
{
PROFILE_START("getVar<unsigned int>");
Array<unsigned int> x( reverse(getVarDim(fid,var)) );
int err = nc_get_var_uint( fid, getVarID(fid,var), x.get() );
CHECK_NC_ERR( err );
PROFILE_STOP("getVar<unsigned int>");
return x;
}
template<>
Array<int> getVar<int>( int fid, const std::string& var )
{
PROFILE_START("getVar<int>");
Array<int> x( reverse(getVarDim(fid,var)) );
int err = nc_get_var_int( fid, getVarID(fid,var), x.get() );
CHECK_NC_ERR( err );
PROFILE_STOP("getVar<int>");
return x;
}
template<>
Array<float> getVar<float>( int fid, const std::string& var )
{
PROFILE_START("getVar<float>");
Array<float> x( reverse(getVarDim(fid,var)) );
int err = nc_get_var_float( fid, getVarID(fid,var), x.get() );
CHECK_NC_ERR( err );
PROFILE_STOP("getVar<float>");
return x;
}
template<>
Array<double> getVar<double>( int fid, const std::string& var )
{
PROFILE_START("getVar<double>");
Array<double> x( reverse(getVarDim(fid,var)) );
int err = nc_get_var_double( fid, getVarID(fid,var), x.get() );
CHECK_NC_ERR( err );
PROFILE_STOP("getVar<double>");
return x;
}
template<>
Array<char> getVar<char>( int fid, const std::string& var )
{
PROFILE_START("getVar<char>");
Array<char> x( reverse(getVarDim(fid,var)) );
int err = nc_get_var_text( fid, getVarID(fid,var), x.get() );
CHECK_NC_ERR( err );
PROFILE_STOP("getVar<char>");
return x;
}
template<>
Array<std::string> getVar<std::string>( int fid, const std::string& var )
{
PROFILE_START("getVar<std::string>");
Array<char> tmp = getVar<char>( fid, var );
std::vector<size_t> dim = tmp.size();
if ( dim.size() == 1 )
dim[0] = 1;
else
dim.erase( dim.begin() );
Array<std::string> text(dim);
for (size_t i=0; i<text.length(); i++)
text(i) = &(tmp(0,i));
PROFILE_STOP("getVar<std::string>");
return text;
}
/****************************************************
* Read an attribute *
****************************************************/
template<>
Array<double> getAtt<double>( int fid, const std::string& att )
{
PROFILE_START("getAtt<double>");
Array<double> x( getAttDim(fid,att) );
int err = nc_get_att_double( fid, NC_GLOBAL, att.c_str(), x.get() );
CHECK_NC_ERR( err );
PROFILE_STOP("getAtt<double>");
return x;
}
template<>
Array<std::string> getAtt<std::string>( int fid, const std::string& att )
{
PROFILE_START("getAtt<std::string>");
char *tmp = new char[getAttDim(fid,att)[0]];
Array<std::string> x(1);
x(0) = tmp;
delete [] tmp;
PROFILE_STOP("getAtt<std::string>");
return x;
}
}; // netcdf namespace
#else
#endif

98
IO/netcdf.h Normal file
View File

@@ -0,0 +1,98 @@
#ifndef NETCDF_READER
#define NETCDF_READER
#include <string>
#include <vector>
#include "common/Array.h"
namespace netcdf {
//! Enum to hold variable type
enum VariableType { BYTE, SHORT, USHORT, INT, UINT, INT64, UINT64, FLOAT, DOUBLE, STRING };
/*!
* @brief Open netcdf file
* @detailed This function opens a netcdf file
* @return This function returns a handle to the file
* @param filename File to open
*/
int open( const std::string& filename );
/*!
* @brief Close netcdf file
* @detailed This function closes a netcdf file
* @param fid Handle to the open file
*/
void close( int fid );
/*!
* @brief Read the variable names
* @detailed This function reads a list of the variable names in the file
* @param fid Handle to the open file
*/
std::vector<std::string> getVarNames( int fid );
/*!
* @brief Read the attribute names
* @detailed This function reads a list of the attribute names in the file
* @param fid Handle to the open file
*/
std::vector<std::string> getAttNames( int fid );
/*!
* @brief Return the variable type
* @detailed This function returns the type for a variable
* @param fid Handle to the open file
* @param var Variable to read
*/
VariableType getVarType( int fid, const std::string& var );
/*!
* @brief Return the attribute type
* @detailed This function returns the type for an attribute
* @param fid Handle to the open file
* @param att Attribute to read
*/
VariableType getAttType( int fid, const std::string& att );
/*!
* @brief Return the variable dimensions
* @detailed This function returns the die for a variable
* @param fid Handle to the open file
* @param var Variable to read
*/
std::vector<size_t> getVarDim( int fid, const std::string& var );
/*!
* @brief Read a variable
* @detailed This function reads a variable with the given name from the file
* @param fid Handle to the open file
* @param var Variable to read
*/
template<class TYPE>
Array<TYPE> getVar( int fid, const std::string& var );
/*!
* @brief Read an attribute
* @detailed This function reads an attribute with the given name from the file
* @param fid Handle to the open file
* @param att Attribute to read
*/
template<class TYPE>
Array<TYPE> getAtt( int fid, const std::string& att );
}; // netcdf namespace
#endif

View File

@@ -75,7 +75,9 @@ ENDMACRO()
# Macro to configure MIC
MACRO( CONFIGURE_MIC )
CHECK_ENABLE_FLAG( USE_MIC 0 )
ADD_DEFINITIONS ( "-D USE_MIC" )
IF ( USE_MIC )
ADD_DEFINITIONS( "-D USE_MIC" )
ENDIF()
ENDMACRO()
@@ -173,6 +175,65 @@ MACRO( CONFIGURE_MPI )
ENDMACRO()
# Macro to find and configure hdf5
MACRO ( CONFIGURE_HDF5 )
CHECK_ENABLE_FLAG( USE_HDF5 0 )
IF ( USE_HDF5 )
# Check if we specified the silo directory
IF ( HDF5_DIRECTORY )
VERIFY_PATH ( ${HDF5_DIRECTORY} )
INCLUDE_DIRECTORIES ( ${HDF5_DIRECTORY}/include )
SET ( HDF5_INCLUDE ${HDF5_DIRECTORY}/include )
FIND_LIBRARY ( HDF5_LIB NAMES hdf5 PATHS ${HDF5_DIRECTORY}/lib NO_DEFAULT_PATH )
FIND_LIBRARY ( HDF5_HL_LIB NAMES hdf5_hl PATHS ${HDF5_DIRECTORY}/lib NO_DEFAULT_PATH )
ELSE()
MESSAGE( FATAL_ERROR "Default search for hdf5 is not yet supported. Use -D HDF5_DIRECTORY=" )
ENDIF()
SET ( HDF5_LIBS
${HDF5_HL_LIB}
${HDF5_LIB}
)
ADD_DEFINITIONS ( "-D USE_HDF5" )
MESSAGE( "Using hdf5" )
MESSAGE( " ${HDF5_LIB}" )
ENDIF()
ENDMACRO()
# Macro to find and configure netcdf
MACRO( CONFIGURE_NETCDF )
CHECK_ENABLE_FLAG( USE_NETCDF 0 )
IF ( USE_NETCDF )
SET( USE_HDF5 1 )
CONFIGURE_HDF5()
IF ( NETCDF_DIRECTORY )
VERIFY_PATH ( ${NETCDF_DIRECTORY} )
INCLUDE_DIRECTORIES ( ${NETCDF_DIRECTORY}/include )
SET ( NETCDF_INCLUDE ${NETCDF_DIRECTORY}/include )
FIND_LIBRARY( NETCDF_NETCDF_LIB NAMES netcdf PATHS ${NETCDF_DIRECTORY}/lib NO_DEFAULT_PATH )
FIND_LIBRARY( NETCDF_HDF5_LIB NAMES hdf5 PATHS ${NETCDF_DIRECTORY}/lib NO_DEFAULT_PATH )
FIND_LIBRARY( NETCDF_HL_LIB NAMES hl PATHS ${NETCDF_DIRECTORY}/lib NO_DEFAULT_PATH )
IF ( NOT NETCDF_NETCDF_LIB )
MESSAGE( FATAL_ERROR "Did not find library for netcdf" )
ENDIF()
SET ( NETCDF_LIBS ${NETCDF_NETCDF_LIB} )
IF ( NETCDF_HDF5_LIB )
SET ( NETCDF_LIBS ${NETCDF_LIBS} ${NETCDF_HDF5_LIB} )
ENDIF()
IF ( NETCDF_HL_LIB )
SET ( NETCDF_LIBS ${NETCDF_LIBS} ${NETCDF_HL_LIB} )
ENDIF()
ELSE()
MESSAGE( FATAL_ERROR "Default search for netcdf is not yet supported. Use -D NETCDF_DIRECTORY=" )
ENDIF()
SET( EXTERNAL_LIBS ${EXTERNAL_LIBS} ${NETCDF_LIBS} ${HDF5_LIBS} )
ADD_DEFINITIONS ( "-D USE_NETCDF" )
MESSAGE( "Using netcdf" )
MESSAGE( " ${NETCDF_LIBS}" )
ENDIF()
ENDMACRO()
# Macro to configure system-specific libraries and flags
MACRO( CONFIGURE_SYSTEM )
# First check/set the compile mode

View File

@@ -1,13 +1,17 @@
# Clear all modules (except modules)
for var in ${LOADEDMODULES//:/ }; do if [ " ${var///*/}" != " modules" ]; then module unload " ${var///*/}" > /dev/null 2>&1; fi; done
#for var in ${LOADEDMODULES//:/ }; do if [ " ${var///*/}" != " modules" ]; then module unload " ${var///*/}" > /dev/null 2>&1; fi; done
source $MODULESHOME/init/bash
module swap PrgEnv-intel PrgEnv-gnu
module load cray-hdf5-parallel
module load cray-netcdf-hdf5parallel
# Set the modules and enviornmental variables
module load DefApps
module load PrgEnv-intel
module unload cray-libsci
module load cray-mpich ugni xpmem udreg altd
module load moab torque
#module load DefApps
#module load PrgEnv-intel
#module unload cray-libsci
#module load cray-mpich ugni xpmem udreg altd
#module load moab torque
module load mercurial git
module load cmake
@@ -29,6 +33,9 @@ cmake \
-D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \
-D USE_CUDA=0 \
-D CUDA_FLAGS="-arch sm_35" \
-D USE_NETCDF=1 \
-D NETCDF_DIRECTORY=$NETCDF_DIR \
-D HDF5_DIRECTORY=$HDF5_DIR \
-D PREFIX=$MEMBERWORK/geo106/eos-LBPM-WIA \
~/LBPM-WIA
#-D CUDA_HOST_COMPILER="/usr/bin/gcc" \

View File

@@ -35,6 +35,9 @@ ADD_LBPM_TEST_PARALLEL( TestMassConservationD3Q7 1 )
ADD_LBPM_TEST_1_2_4( testCommunication )
ADD_LBPM_TEST_1_2_4( testUtilities )
ADD_LBPM_TEST( TestWriter )
IF ( USE_NETCDF )
ADD_LBPM_PROVISIONAL_TEST( TestNetcdf )
ENDIF()
# Sample test that will run with 1, 2, and 4 processors, failing with 4 or more procs
ADD_LBPM_TEST_1_2_4( hello_world )

57
tests/TestNetcdf.cpp Normal file
View File

@@ -0,0 +1,57 @@
// Sequential blob analysis
// Reads parallel simulation data and performs connectivity analysis
// and averaging on a blob-by-blob basis
// James E. McClure 2014
#include "IO/netcdf.h"
#include "ProfilerApp.h"
void load( const std::string filename )
{
int fid = netcdf::open( filename );
std::vector<std::string> vars = netcdf::getVarNames( fid );
for (size_t i=0; i<vars.size(); i++) {
printf("Reading variable %s\n",vars[i].c_str());
netcdf::VariableType type = netcdf::getVarType( fid, vars[i] );
if ( type == netcdf::STRING )
Array<std::string> tmp = netcdf::getVar<std::string>( fid, vars[i] );
else if ( type == netcdf::SHORT )
Array<short> tmp = netcdf::getVar<short>( fid, vars[i] );
else
Array<double> tmp = netcdf::getVar<double>( fid, vars[i] );
}
std::vector<std::string> attr = netcdf::getAttNames( fid );
for (size_t i=0; i<attr.size(); i++) {
printf("Reading attribute %s\n",attr[i].c_str());
netcdf::VariableType type = netcdf::getAttType( fid, attr[i] );
if ( type == netcdf::STRING )
Array<std::string> tmp = netcdf::getAtt<std::string>( fid, attr[i] );
else
Array<double> tmp = netcdf::getAtt<double>( fid, attr[i] );
}
netcdf::close( fid );
}
int main(int argc, char **argv)
{
PROFILE_START("Main");
std::vector<std::string> filenames;
if ( argc==0 ) {
printf("At least one filename must be specified\n");
return 1;
}
for (int i=1; i<argc; i++)
load( argv[i] );
PROFILE_SAVE("TestNetcdf");
return 0;
}

View File

@@ -21,7 +21,14 @@ int main(int argc, char **argv)
MPI_Comm comm = MPI_COMM_WORLD;
MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs);
int SOLID=atoi(argv[1]);
int NWP=atoi(argv[2]);
if (rank==0){
printf("Solid Label %i \n",SOLID);
printf("NWP Label %i \n",NWP);
}
//.......................................................................
// Reading the domain information file
//.......................................................................
@@ -183,14 +190,18 @@ int main(int argc, char **argv)
int count = 0;
N=nx*ny*nz;
if (rank==0) printf("WARNING: assumed that INPUT: WP(1),NWP(2) OUTPUT will be NWP(1),WP(2) (lbpm_segmented_decomp) \n");
//if (rank==0) printf("WARNING: assumed that INPUT: WP(1),NWP(2) OUTPUT will be NWP(1),WP(2) (lbpm_segmented_decomp) \n");
// Really need a better way to do this -- this is to flip convention for a particular data set
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;
if (Dm.id[n] == 1) Dm.id[n]=2;
else if (Dm.id[n] == 2) Dm.id[n]=1;
n = k*nx*ny+j*nx+i;
if (Dm.id[n]==char(SOLID)) Dm.id[n] = 0;
else if (Dm.id[n]==char(NWP)) Dm.id[n] = 1;
else Dm.id[n] = 2;
// if (Dm.id[n] == 1) Dm.id[n]=2;
//else if (Dm.id[n] == 2) Dm.id[n]=1;
// Initialize the solid phase
// if (Dm.id[n] == 0) id[n] = 0;
// else id[n] = 1;
@@ -198,6 +209,27 @@ int main(int argc, char **argv)
}
}
count = 0;
int total = 0;
int countGlobal = 0;
int totalGlobal = 0;
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
n=k*Nx*Ny+j*Nx+i;
total++;
if (Dm.id[n] == 0){
count++;
}
}
}
}
MPI_Allreduce(&count,&countGlobal,1,MPI_INT,MPI_SUM,comm);
MPI_Allreduce(&total,&totalGlobal,1,MPI_INT,MPI_SUM,comm);
float porosity = float(totalGlobal-countGlobal)/totalGlobal;
if (rank==0) printf("Porosity=%f\n",porosity);
char LocalRankFilename[40];
sprintf(LocalRankFilename,"ID.%05i",rank);