From b5da2246249ade3b6472b6bbbad1e6c018a8685f Mon Sep 17 00:00:00 2001
From: Mark Berrill <berrillma@ornl.gov>
Date: Mon, 22 Feb 2016 16:26:35 -0500
Subject: [PATCH] Creating initial reader for netcdf

---
 CMakeLists.txt        |   1 +
 IO/netcdf.cpp         | 288 ++++++++++++++++++++++++++++++++++++++++++
 IO/netcdf.h           |  98 ++++++++++++++
 cmake/libraries.cmake |  63 ++++++++-
 tests/CMakeLists.txt  |   3 +
 tests/TestNetcdf.cpp  |  57 +++++++++
 6 files changed, 509 insertions(+), 1 deletion(-)
 create mode 100644 IO/netcdf.cpp
 create mode 100644 IO/netcdf.h
 create mode 100644 tests/TestNetcdf.cpp

diff --git a/CMakeLists.txt b/CMakeLists.txt
index b67cf8ca..011f5ab4 100755
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -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()
diff --git a/IO/netcdf.cpp b/IO/netcdf.cpp
new file mode 100644
index 00000000..6c52b049
--- /dev/null
+++ b/IO/netcdf.cpp
@@ -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
+
+
diff --git a/IO/netcdf.h b/IO/netcdf.h
new file mode 100644
index 00000000..5eda1037
--- /dev/null
+++ b/IO/netcdf.h
@@ -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
diff --git a/cmake/libraries.cmake b/cmake/libraries.cmake
index 9ed75e5a..bb32f16b 100644
--- a/cmake/libraries.cmake
+++ b/cmake/libraries.cmake
@@ -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
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index d1e20511..50a374df 100755
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -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 )
diff --git a/tests/TestNetcdf.cpp b/tests/TestNetcdf.cpp
new file mode 100644
index 00000000..cf8e55aa
--- /dev/null
+++ b/tests/TestNetcdf.cpp
@@ -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;  
+}
+