From 90ffb5b9d062287d3f01cfcf6e55b6bd862ba3b1 Mon Sep 17 00:00:00 2001
From: Mark Berrill <berrillma@ornl.gov>
Date: Wed, 1 Mar 2017 15:48:17 -0500
Subject: [PATCH] Adding support for different precision in silo

---
 IO/Mesh.cpp                          |  40 +++---
 IO/Mesh.h                            |  20 +--
 IO/Reader.cpp                        |  24 ++--
 IO/Writer.cpp                        |  87 ++++++++++--
 IO/silo.cpp                          | 191 ++++++++++++++++++---------
 IO/silo.h                            |  45 +++++--
 tests/ComponentLabel.cpp             |  10 +-
 tests/TestBlobIdentify.cpp           |   6 +-
 tests/TestBubble.cpp                 |   2 +-
 tests/TestWriter.cpp                 |  81 +++++++-----
 tests/lbpm_color_macro_simulator.cpp |   8 +-
 tests/lbpm_color_simulator.cpp       |   8 +-
 tests/lbpm_uCT_maskfilter.cpp        |  16 +--
 tests/lbpm_uCT_pp.cpp                |  16 +--
 14 files changed, 366 insertions(+), 188 deletions(-)

diff --git a/IO/Mesh.cpp b/IO/Mesh.cpp
index f4336117..742dac85 100644
--- a/IO/Mesh.cpp
+++ b/IO/Mesh.cpp
@@ -48,13 +48,13 @@ bool MeshDataStruct::check() const
         if ( mesh2 == nullptr )
             return false;
         for ( const auto& var : vars ) {
-            if ( var->type == IO::NodeVariable ) {
+            if ( var->type == IO::VariableType::NodeVariable ) {
                 pass = pass && var->data.size(0)==mesh2->points.size() && var->data.size(1)==var->dim;
-            } else if ( var->type == IO::EdgeVariable ) {
+            } else if ( var->type == IO::VariableType::EdgeVariable ) {
                 ERROR("Invalid type for PointList");
-            } else if ( var->type == IO::SurfaceVariable ) {
+            } else if ( var->type == IO::VariableType::SurfaceVariable ) {
                 ERROR("Invalid type for PointList");
-            } else if ( var->type == IO::VolumeVariable ) {
+            } else if ( var->type == IO::VariableType::VolumeVariable ) {
                 ERROR("Invalid type for PointList");
             } else {
                 ERROR("Invalid variable type");
@@ -65,13 +65,13 @@ bool MeshDataStruct::check() const
         if ( mesh2 == nullptr )
             return false;
         for ( const auto& var : vars ) {
-            if ( var->type == IO::NodeVariable ) {
+            if ( var->type == IO::VariableType::NodeVariable ) {
                 pass = pass && var->data.size(0)==mesh2->vertices->points.size() && var->data.size(1)==var->dim;
-            } else if ( var->type == IO::EdgeVariable ) {
+            } else if ( var->type == IO::VariableType::EdgeVariable ) {
                 ERROR("Not finished");
-            } else if ( var->type == IO::SurfaceVariable ) {
+            } else if ( var->type == IO::VariableType::SurfaceVariable ) {
                 ERROR("Not finished");
-            } else if ( var->type == IO::VolumeVariable ) {
+            } else if ( var->type == IO::VariableType::VolumeVariable ) {
                 pass = pass && var->data.size(0)==mesh2->A.size() && var->data.size(1)==var->dim;
             } else {
                 ERROR("Invalid variable type");
@@ -82,14 +82,14 @@ bool MeshDataStruct::check() const
         if ( mesh2 == nullptr )
             return false;
         for ( const auto& var : vars ) {
-            if ( var->type == IO::NodeVariable ) {
+            if ( var->type == IO::VariableType::NodeVariable ) {
                 pass = pass && (int) var->data.size(0)==(mesh2->nx+1) && (int) var->data.size(1)==(mesh2->ny+1)
                             && (int) var->data.size(2)==(mesh2->nz+1) &&       var->data.size(3)==var->dim;
-            } else if ( var->type == IO::EdgeVariable ) {
+            } else if ( var->type == IO::VariableType::EdgeVariable ) {
                 ERROR("Not finished");
-            } else if ( var->type == IO::SurfaceVariable ) {
+            } else if ( var->type == IO::VariableType::SurfaceVariable ) {
                 ERROR("Not finished");
-            } else if ( var->type == IO::VolumeVariable ) {
+            } else if ( var->type == IO::VariableType::VolumeVariable ) {
                 pass = pass && (int) var->data.size(0)==mesh2->nx && (int) var->data.size(1)==mesh2->ny
                             && (int) var->data.size(2)==mesh2->nz &&       var->data.size(3)==var->dim;
             } else {
@@ -120,7 +120,7 @@ PointList::~PointList( )
 size_t PointList::numberPointsVar( VariableType type ) const
 {
     size_t N = 0;
-    if ( type == NodeVariable )
+    if ( type == VariableType::NodeVariable )
         N = points.size();
     return N;
 }
@@ -196,9 +196,9 @@ TriList::~TriList( )
 size_t TriList::numberPointsVar( VariableType type ) const
 {
     size_t N = 0;
-    if ( type==NodeVariable )
+    if ( type==VariableType::NodeVariable )
         N = 3*A.size();
-    else if ( type==SurfaceVariable ||  type==VolumeVariable )
+    else if ( type==VariableType::SurfaceVariable ||  type==VariableType::VolumeVariable )
         N = A.size();
     return N;
 }
@@ -301,9 +301,9 @@ TriMesh::~TriMesh( )
 size_t TriMesh::numberPointsVar( VariableType type ) const
 {
     size_t N = 0;
-    if ( type==NodeVariable )
+    if ( type==VariableType::NodeVariable )
         N = vertices->points.size();
-    else if ( type==SurfaceVariable ||  type==VolumeVariable )
+    else if ( type==VariableType::SurfaceVariable ||  type==VariableType::VolumeVariable )
         N = A.size();
     return N;
 }
@@ -387,11 +387,11 @@ DomainMesh::~DomainMesh()
 size_t DomainMesh::numberPointsVar( VariableType type ) const
 {
     size_t N = 0;
-    if ( type==NodeVariable )
+    if ( type==VariableType::NodeVariable )
         N = (nx+1)*(ny+1)*(nz+1);
-    else if ( type==SurfaceVariable )
+    else if ( type==VariableType::SurfaceVariable )
         N = (nx+1)*ny*nz + nx*(ny+1)*nz + nx*ny*(nz+1);
-    else if ( type==VolumeVariable )
+    else if ( type==VariableType::VolumeVariable )
         N = nx*ny*nz;
     return N;
 }
diff --git a/IO/Mesh.h b/IO/Mesh.h
index f2ff88b7..969a9a61 100644
--- a/IO/Mesh.h
+++ b/IO/Mesh.h
@@ -16,8 +16,8 @@ namespace IO {
 
 
 //! Possible variable types
-//enum class VariableType : unsigned char { NodeVariable=1, EdgeVariable=2, SurfaceVariable=2, VolumeVariable=3, Null=0 };
-enum VariableType { NodeVariable=1, EdgeVariable=2, SurfaceVariable=2, VolumeVariable=3, NullVariable=0 };
+enum class VariableType: unsigned char { NodeVariable=1, EdgeVariable=2, SurfaceVariable=2, VolumeVariable=3, NullVariable=0 };
+enum class DataType: unsigned char { Double=1, Float=2, Int=2, Null=0 };
 
 
 /*! \class Mesh
@@ -168,18 +168,19 @@ struct Variable
 {
 public:
     // Internal variables
-    unsigned int dim;           //!< Number of points per grid point (1: scalar, 3: vector, ...)
+    unsigned char dim;          //!< Number of points per grid point (1: scalar, 3: vector, ...)
     VariableType type;          //!< Variable type
+    DataType precision;         //!< Variable precision to use for IO
     std::string name;           //!< Variable name
     Array<double> data;         //!< Variable data
     //! Empty constructor
-    Variable(): dim(0), type(NullVariable) {}
+    Variable(): dim(0), type(VariableType::NullVariable), precision(DataType::Double) {}
     //! Constructor
     Variable( int dim_, IO::VariableType type_, const std::string& name_ ):
-        dim(dim_), type(type_), name(name_) {}
+        dim(dim_), type(type_), precision(DataType::Double), name(name_) {}
     //! Constructor
     Variable( int dim_, IO::VariableType type_, const std::string& name_, const Array<double>& data_ ):
-        dim(dim_), type(type_), name(name_), data(data_) {}
+        dim(dim_), type(type_), precision(DataType::Double), name(name_), data(data_) {}
     //! Destructor
     virtual ~Variable() {}
 protected:
@@ -194,9 +195,12 @@ protected:
     \brief A class used to hold database info for saving a mesh
 */
 struct MeshDataStruct {
-    std::string             meshName;
-    std::shared_ptr<Mesh>   mesh;
+    DataType precision;         //!< Precision to use for IO (mesh)
+    std::string meshName;       //!< Mesh name
+    std::shared_ptr<Mesh> mesh; //!< Mesh data
     std::vector<std::shared_ptr<Variable> >  vars;
+    //! Empty constructor
+    MeshDataStruct(): precision(DataType::Double) {}
     //! Check the data
     bool check() const;
 };
diff --git a/IO/Reader.cpp b/IO/Reader.cpp
index cf60c388..c9595a8b 100644
--- a/IO/Reader.cpp
+++ b/IO/Reader.cpp
@@ -163,7 +163,7 @@ std::shared_ptr<IO::Mesh> IO::getMesh( const std::string& path, const std::strin
         int rank = std::stoi(database.file.substr(0,database.file.find(".silo")).c_str());
         auto fid = silo::open( filename, silo::READ );
         if ( meshDatabase.meshClass=="PointList" ) {
-            Array<double> coords = silo::readPointMesh( fid, database.name );
+            Array<double> coords = silo::readPointMesh<double>( fid, database.name );
             ASSERT(coords.size(1)==3);
             std::shared_ptr<IO::PointList> mesh2( new IO::PointList( coords.size(0) ) );
             for (size_t i=0; i<coords.size(1); i++) {
@@ -266,11 +266,11 @@ std::shared_ptr<IO::Variable> IO::getVariable( const std::string& path, const st
         auto fid = silo::open( filename, silo::READ );
         var.reset( new Variable( variableDatabase.dim, variableDatabase.type, variable ) );
         if ( meshDatabase.meshClass=="PointList" ) {
-            var->data = silo::readPointMeshVariable( fid, variable );
+            var->data = silo::readPointMeshVariable<double>( fid, variable );
         } else if ( meshDatabase.meshClass=="TriMesh" || meshDatabase.meshClass=="TriList" ) {
-            var->data = silo::readTriMeshVariable( fid, variable );
+            var->data = silo::readTriMeshVariable<double>( fid, variable );
         } else if ( meshDatabase.meshClass=="DomainMesh" ) {
-            var->data = silo::readUniformMeshVariable( fid, variable );
+            var->data = silo::readUniformMeshVariable<double>( fid, variable );
         } else {
             ERROR("Unknown mesh class");
         }
@@ -293,15 +293,15 @@ void IO::reformatVariable( const IO::Mesh& mesh, IO::Variable& var )
 {
     if ( mesh.className() == "DomainMesh" ) {
         const IO::DomainMesh& mesh2 = dynamic_cast<const IO::DomainMesh&>( mesh );
-        if ( var.type == NodeVariable ) {
+        if ( var.type == VariableType::NodeVariable ) {
             size_t N2 = var.data.length() / ((mesh2.nx+1)*(mesh2.ny+1)*(mesh2.nz+1));
             ASSERT( (mesh2.nx+1)*(mesh2.ny+1)*(mesh2.nz+1)*N2 == var.data.length() );
             var.data.reshape( { (size_t) mesh2.nx+1, (size_t) mesh2.ny+1, (size_t) mesh2.nz+1, N2 } );
-        } else if ( var.type == EdgeVariable ) {
+        } else if ( var.type == VariableType::EdgeVariable ) {
             ERROR("Not finished");
-        } else if ( var.type == SurfaceVariable ) {
+        } else if ( var.type == VariableType::SurfaceVariable ) {
             ERROR("Not finished");
-        } else if ( var.type == VolumeVariable ) {
+        } else if ( var.type == VariableType::VolumeVariable ) {
             size_t N2 = var.data.length() / (mesh2.nx*mesh2.ny*mesh2.nz);
             ASSERT( mesh2.nx*mesh2.ny*mesh2.nz*N2 == var.data.length() );
             var.data.reshape( { (size_t) mesh2.nx, (size_t) mesh2.ny, (size_t) mesh2.nz, N2 } );
@@ -317,16 +317,16 @@ void IO::reformatVariable( const IO::Mesh& mesh, IO::Variable& var )
     } else if ( mesh.className()=="TriMesh" || mesh.className() == "TriList" ) {
         std::shared_ptr<Mesh> mesh_ptr( const_cast<Mesh*>(&mesh), []( void* ) {} );
         std::shared_ptr<TriMesh> mesh2 = getTriMesh( mesh_ptr );
-        if ( var.type == NodeVariable ) {
+        if ( var.type == VariableType::NodeVariable ) {
             size_t N = mesh2->vertices->points.size();
             size_t N_var = var.data.length()/N;
             ASSERT( N*N_var == var.data.length() );
             var.data.reshape( { N, N_var } );
-        } else if ( var.type == EdgeVariable ) {
+        } else if ( var.type == VariableType::EdgeVariable ) {
             ERROR("Not finished");
-        } else if ( var.type == SurfaceVariable ) {
+        } else if ( var.type == VariableType::SurfaceVariable ) {
             ERROR("Not finished");
-        } else if ( var.type == VolumeVariable ) {
+        } else if ( var.type == VariableType::VolumeVariable ) {
             size_t N = mesh2->A.size();
             size_t N_var = var.data.length()/N;
             ASSERT( N*N_var == var.data.length() );
diff --git a/IO/Writer.cpp b/IO/Writer.cpp
index 10ef6b24..12bf07d0 100644
--- a/IO/Writer.cpp
+++ b/IO/Writer.cpp
@@ -180,7 +180,7 @@ static IO::MeshDatabase write_domain( FILE *fid, const std::string& filename,
         int dim = mesh.vars[i]->dim;
         int type = static_cast<int>(mesh.vars[i]->type);
         size_t N = mesh.vars[i]->data.length();
-        if ( type == static_cast<int>(IO::NullVariable) ) {
+        if ( type == static_cast<int>(IO::VariableType::NullVariable) ) {
             ERROR("Variable type not set");
         }
         size_t N_mesh = mesh.mesh->numberPointsVar(mesh.vars[i]->type);
@@ -197,10 +197,30 @@ static IO::MeshDatabase write_domain( FILE *fid, const std::string& filename,
 
 #ifdef USE_SILO
 // Write a PointList mesh (and variables) to a file
+template<class TYPE>
+static void writeSiloPointMesh( DBfile *fid, const IO::PointList& mesh, const std::string& meshname )
+{
+    const auto& points = mesh.getPoints();
+    std::vector<TYPE> x(points.size()), y(points.size()), z(points.size());
+    for (size_t i=0; i<x.size(); i++) {
+        x[i] = points[i].x;
+        y[i] = points[i].y;
+        z[i] = points[i].z;
+    }
+    const TYPE *coords[] = { x.data(), y.data(), z.data() };
+    silo::writePointMesh<TYPE>( fid, meshname, 3, points.size(), coords );
+}
 static void writeSiloPointList( DBfile *fid, const IO::MeshDataStruct& meshData, IO::MeshDatabase database )
 {
     const IO::PointList& mesh = dynamic_cast<IO::PointList&>( *meshData.mesh );
     const std::string meshname = database.domains[0].name;
+    if ( meshData.precision == IO::DataType::Double ) {
+        writeSiloPointMesh<double>( fid, mesh, meshname );
+    } else if ( meshData.precision == IO::DataType::Float ) {
+        writeSiloPointMesh<float>( fid, mesh, meshname );
+    } else {
+        ERROR("Unsupported format");
+    }
     const auto& points = mesh.getPoints();
     std::vector<double> x(points.size()), y(points.size()), z(points.size());
     for (size_t i=0; i<x.size(); i++) {
@@ -212,28 +232,63 @@ static void writeSiloPointList( DBfile *fid, const IO::MeshDataStruct& meshData,
     silo::writePointMesh( fid, meshname, 3, points.size(), coords );
     for (size_t i=0; i<meshData.vars.size(); i++) {
         const IO::Variable& var = *meshData.vars[i];
-        silo::writePointMeshVariable( fid, meshname, var.name, var.data );
+        if ( var.precision == IO::DataType::Double ) {
+            silo::writePointMeshVariable( fid, meshname, var.name, var.data );
+        } else if ( var.precision == IO::DataType::Float ) {
+            Array<float> data2( var.data.size() );
+            data2.copy( var.data );
+            silo::writePointMeshVariable( fid, meshname, var.name, data2 );
+        } else if ( var.precision == IO::DataType::Int ) {
+            Array<int> data2( var.data.size() );
+            data2.copy( var.data );
+            silo::writePointMeshVariable( fid, meshname, var.name, data2 );
+        } else {
+            ERROR("Unsupported format");
+        }
     }
 }
 // Write a TriMesh mesh (and variables) to a file
-static void writeSiloTriMesh2( DBfile *fid, const IO::MeshDataStruct& meshData,
-    const IO::TriMesh& mesh, IO::MeshDatabase database )
+template<class TYPE>
+static void writeSiloTriMesh( DBfile *fid, const IO::TriMesh& mesh, const std::string& meshname )
 {
-    const std::string meshname = database.domains[0].name;
     const auto& points = mesh.vertices->getPoints();
-    std::vector<double> x(points.size()), y(points.size()), z(points.size());
+    std::vector<TYPE> x(points.size()), y(points.size()), z(points.size());
     for (size_t i=0; i<x.size(); i++) {
         x[i] = points[i].x;
         y[i] = points[i].y;
         z[i] = points[i].z;
     }
-    const double *coords[] = { x.data(), y.data(), z.data() };
+    const TYPE *coords[] = { x.data(), y.data(), z.data() };
     const int *tri[] = { mesh.A.data(), mesh.B.data(), mesh.C.data() };
-    silo::writeTriMesh( fid, meshname, 3, 2, points.size(), coords, mesh.A.size(), tri );
+    silo::writeTriMesh<TYPE>( fid, meshname, 3, 2, points.size(), coords, mesh.A.size(), tri );
+}
+static void writeSiloTriMesh2( DBfile *fid, const IO::MeshDataStruct& meshData,
+    const IO::TriMesh& mesh, IO::MeshDatabase database )
+{
+    const std::string meshname = database.domains[0].name;
+    if ( meshData.precision == IO::DataType::Double ) {
+        writeSiloTriMesh<double>( fid, mesh, meshname );
+    } else if ( meshData.precision == IO::DataType::Float ) {
+        writeSiloTriMesh<float>( fid, mesh, meshname );
+    } else {
+        ERROR("Unsupported format");
+    }
     for (size_t i=0; i<meshData.vars.size(); i++) {
         const IO::Variable& var = *meshData.vars[i];
         auto type = static_cast<silo::VariableType>( var.type );
-        silo::writeTriMeshVariable( fid, 3, meshname, var.name, var.data, type );
+        if ( var.precision == IO::DataType::Double ) {
+            silo::writeTriMeshVariable( fid, 3, meshname, var.name, var.data, type );
+        } else if ( var.precision == IO::DataType::Float ) {
+            Array<float> data2( var.data.size() );
+            data2.copy( var.data );
+            silo::writeTriMeshVariable( fid, 3, meshname, var.name, data2, type );
+        } else if ( var.precision == IO::DataType::Int ) {
+            Array<int> data2( var.data.size() );
+            data2.copy( var.data );
+            silo::writeTriMeshVariable( fid, 3, meshname, var.name, data2, type );
+        } else {
+            ERROR("Unsupported format");
+        }
     }
 }
 static void writeSiloTriMesh( DBfile *fid, const IO::MeshDataStruct& meshData, IO::MeshDatabase database )
@@ -262,7 +317,19 @@ static void writeSiloDomainMesh( DBfile *fid, const IO::MeshDataStruct& meshData
     for (size_t i=0; i<meshData.vars.size(); i++) {
         const IO::Variable& var = *meshData.vars[i];
         auto type = static_cast<silo::VariableType>( var.type );
-        silo::writeUniformMeshVariable<3>( fid, meshname, N, var.name, var.data, type );
+        if ( var.precision == IO::DataType::Double ) {
+            silo::writeUniformMeshVariable<3>( fid, meshname, N, var.name, var.data, type );
+        } else if ( var.precision == IO::DataType::Float ) {
+            Array<float> data2( var.data.size() );
+            data2.copy( var.data );
+            silo::writeUniformMeshVariable<3>( fid, meshname, N, var.name, data2, type );
+        } else if ( var.precision == IO::DataType::Int ) {
+            Array<int> data2( var.data.size() );
+            data2.copy( var.data );
+            silo::writeUniformMeshVariable<3>( fid, meshname, N, var.name, data2, type );
+        } else {
+            ERROR("Unsupported format");
+        }
     }
 }
 // Write a mesh (and variables) to a file
diff --git a/IO/silo.cpp b/IO/silo.cpp
index 03fde877..dfcd626c 100644
--- a/IO/silo.cpp
+++ b/IO/silo.cpp
@@ -14,6 +14,41 @@
 namespace silo {
 
 
+/****************************************************
+* Helper functions                                  *
+****************************************************/
+template<class TYPE> static constexpr int getType();
+template<> constexpr int getType<double>() { return DB_DOUBLE; }
+template<> constexpr int getType<float>() { return DB_FLOAT; }
+template<> constexpr int getType<int>() { return DB_INT; }
+VariableDataType varDataType( DBfile *fid, const std::string& name )
+{
+    auto type = DBGetVarType( fid, name.c_str() );
+    VariableDataType type2 = VariableDataType::UNKNOWN;
+    if ( type == DB_DOUBLE )
+        type2 = VariableDataType::DOUBLE;
+    else if ( type == DB_FLOAT )
+        type2 = VariableDataType::FLOAT;
+    else if ( type == DB_INT )
+        type2 = VariableDataType::INT;
+    return type2;
+}
+template<class TYPE>
+static inline void copyData( Array<TYPE>& data, int type, const void *src )
+{
+    if ( type == getType<TYPE>() )
+        memcpy( data.data(), src, data.length()*sizeof(TYPE) );
+    else if ( type == DB_DOUBLE )
+        data.copy<double>( static_cast<const double*>(src) );
+    else if ( type == DB_FLOAT )
+        data.copy<float>( static_cast<const float*>(src) );
+    else if ( type == DB_INT )
+        data.copy<int>( static_cast<const int*>(src) );
+    else
+        ERROR("Unknown type");
+}
+
+
 /****************************************************
 * Open/close a file                                 *
 ****************************************************/
@@ -162,15 +197,15 @@ void readUniformMesh( DBfile* fid, const std::string& meshname,
 /****************************************************
 * Write a vector/tensor quad variable               *
 ****************************************************/
-template<int NDIM>
+template<int NDIM,class TYPE>
 void writeUniformMeshVariable( DBfile* fid, const std::string& meshname, const std::array<int,NDIM>& N,
-    const std::string& varname, const Array<double>& data, VariableType type )
+    const std::string& varname, const Array<TYPE>& data, VariableType type )
 {
     PROFILE_START("writeUniformMeshVariable",2);
     int nvars=1, dims[NDIM]={1};
-    const double *vars[NDIM] = { nullptr };
+    const TYPE *vars[NDIM] = { nullptr };
     int vartype = 0;
-    if ( type == NodeVariable ) {
+    if ( type == VariableType::NodeVariable ) {
         ASSERT( data.ndim()==NDIM || data.ndim()==NDIM+1 );
         for (int d=0; d<NDIM; d++)
             ASSERT(N[d]+1==(int)data.size(d));
@@ -181,11 +216,11 @@ void writeUniformMeshVariable( DBfile* fid, const std::string& meshname, const s
             dims[d] = data.size(d);
         for (int i=0; i<nvars; i++)
             vars[i] = &data(i*N);
-    } else if ( type == EdgeVariable ) {
+    } else if ( type == VariableType::EdgeVariable ) {
         ERROR("Not finished");
-    } else if ( type == SurfaceVariable ) {
+    } else if ( type == VariableType::SurfaceVariable ) {
         ERROR("Not finished");
-    } else if ( type == VolumeVariable ) {
+    } else if ( type == VariableType::VolumeVariable ) {
         ASSERT( data.ndim()==NDIM || data.ndim()==NDIM+1 );
         for (int d=0; d<NDIM; d++)
             ASSERT(N[d]==(int)data.size(d));
@@ -207,20 +242,21 @@ void writeUniformMeshVariable( DBfile* fid, const std::string& meshname, const s
     for (int i=0; i<nvars; i++)
         varnames[i] = const_cast<char*>(var_names[i].c_str());
     int err = DBPutQuadvar( fid, varname.c_str(), meshname.c_str(), nvars,
-        varnames.data(), vars, dims, NDIM, nullptr, 0, DB_DOUBLE, vartype, nullptr );
+        varnames.data(), vars, dims, NDIM, nullptr, 0, getType<TYPE>(), vartype, nullptr );
     ASSERT( err == 0 );
     PROFILE_STOP("writeUniformMeshVariable",2);
 }
-Array<double> readUniformMeshVariable( DBfile* fid, const std::string& varname )
+template<class TYPE> 
+Array<TYPE> readUniformMeshVariable( DBfile* fid, const std::string& varname )
 {
     auto var = DBGetQuadvar( fid, varname.c_str() );
     ASSERT( var != nullptr );
-    Array<double> data( var->nels, var->nvals );
-    if ( var->datatype == DB_DOUBLE ) {
-        for (int i=0; i<var->nvals; i++)
-            memcpy( &data(0,i), var->vals[i], var->nels*sizeof(double) );
-    } else {
-        ERROR("Unsupported format");
+    Array<TYPE> data( var->nels, var->nvals );
+    int type = var->datatype;
+    for (int i=0; i<var->nvals; i++) {
+        Array<TYPE> data2( var->nels );
+        copyData<TYPE>( data2, type, var->vals[i] );
+        memcpy( &data(0,i), data2.data(), var->nels*sizeof(TYPE) );
     }
     DBFreeQuadvar( var );
     std::vector<size_t> dims( var->ndims+1, var->nvals );
@@ -234,48 +270,52 @@ Array<double> readUniformMeshVariable( DBfile* fid, const std::string& varname )
 /****************************************************
 * Read/write a point mesh/variable to silo          *
 ****************************************************/
+template<class TYPE>
 void writePointMesh( DBfile* fid, const std::string& meshname,
-    int ndim, int N, const double *coords[] )
+    int ndim, int N, const TYPE *coords[] )
 {
-    int err = DBPutPointmesh( fid, meshname.c_str(), ndim, coords, N, DB_DOUBLE, nullptr );
+    int err = DBPutPointmesh( fid, meshname.c_str(), ndim, coords, N, getType<TYPE>(), nullptr );
     ASSERT( err == 0 );
 }
-Array<double> readPointMesh( DBfile* fid, const std::string& meshname )
+template<class TYPE> 
+Array<TYPE> readPointMesh( DBfile* fid, const std::string& meshname )
 {
     auto mesh = DBGetPointmesh( fid, meshname.c_str() );
     int N = mesh->nels;
     int ndim = mesh->ndims;
-    Array<double> coords(N,ndim);
-    if ( mesh->datatype == DB_DOUBLE ) {
-        for (int d=0; d<ndim; d++)
-            memcpy(&coords(0,d),mesh->coords[d],N*sizeof(double));
-    } else {
-        ERROR("Unsupported format");
+    Array<TYPE> coords(N,ndim);
+    int type = mesh->datatype;
+    for (int d=0; d<ndim; d++) {
+        Array<TYPE> data2( N );
+        copyData<TYPE>( data2, type, mesh->coords[d] );
+        memcpy( &coords(0,d), data2.data(), N*sizeof(TYPE) );
     }
     DBFreePointmesh( mesh );
     return coords;
 }
+template<class TYPE>
 void writePointMeshVariable( DBfile* fid, const std::string& meshname,
-    const std::string& varname, const Array<double>& data )
+    const std::string& varname, const Array<TYPE>& data )
 {
     int N = data.size(0);
     int nvars = data.size(1);
-    std::vector<const double*> vars(nvars);
+    std::vector<const TYPE*> vars(nvars);
     for (int i=0; i<nvars; i++)
         vars[i] = &data(0,i);
-    int err = DBPutPointvar( fid, varname.c_str(), meshname.c_str(), nvars, vars.data(), N, DB_DOUBLE, nullptr );
+    int err = DBPutPointvar( fid, varname.c_str(), meshname.c_str(), nvars, vars.data(), N, getType<TYPE>(), nullptr );
     ASSERT( err == 0 );
 }
-Array<double> readPointMeshVariable( DBfile* fid, const std::string& varname )
+template<class TYPE> 
+Array<TYPE> readPointMeshVariable( DBfile* fid, const std::string& varname )
 {
     auto var = DBGetPointvar( fid, varname.c_str() );
     ASSERT( var != nullptr );
-    Array<double> data( var->nels, var->nvals );
-    if ( var->datatype == DB_DOUBLE ) {
-        for (int i=0; i<var->nvals; i++)
-            memcpy( &data(0,i), var->vals[i], var->nels*sizeof(double) );
-    } else {
-        ERROR("Unsupported format");
+    Array<TYPE> data( var->nels, var->nvals );
+    int type = var->datatype;
+    for (int i=0; i<var->nvals; i++) {
+        Array<TYPE> data2( var->nels );
+        copyData<TYPE>( data2, type, var->vals[i] );
+        memcpy( &data(0,i), data2.data(), var->nels*sizeof(TYPE) );
     }
     DBFreeMeshvar( var );
     return data;
@@ -285,8 +325,9 @@ Array<double> readPointMeshVariable( DBfile* fid, const std::string& varname )
 /****************************************************
 * Read/write a triangle mesh                        *
 ****************************************************/
+template<class TYPE>
 void writeTriMesh( DBfile* fid, const std::string& meshName,
-    int ndim, int ndim_tri, int N, const double *coords[], int N_tri, const int *tri[] )
+    int ndim, int ndim_tri, int N, const TYPE *coords[], int N_tri, const int *tri[] )
 {
     auto zoneName = meshName + "_zones";
     std::vector<int> nodelist( (ndim_tri+1)*N_tri );
@@ -308,25 +349,26 @@ void writeTriMesh( DBfile* fid, const std::string& meshName,
     DBPutZonelist2( fid, zoneName.c_str(), N_tri, ndim_tri, nodelist.data(),
         nodelist.size(), 0, 0, 0, &shapetype, &shapesize, &shapecnt, 1, nullptr );
     DBPutUcdmesh( fid, meshName.c_str(), ndim, nullptr, coords, N,
-        nodelist.size(), zoneName.c_str(), nullptr, DB_DOUBLE, nullptr );
+        nodelist.size(), zoneName.c_str(), nullptr, getType<TYPE>(), nullptr );
 }
-void readTriMesh( DBfile* fid, const std::string& meshname, Array<double>& coords, Array<int>& tri )
+template<class TYPE>
+void readTriMesh( DBfile* fid, const std::string& meshname, Array<TYPE>& coords, Array<int>& tri )
 {
     auto mesh = DBGetUcdmesh( fid, meshname.c_str() );
     int ndim = mesh->ndims;
     int N_nodes = mesh->nnodes;
     coords.resize(N_nodes,ndim);
-    if ( mesh->datatype == DB_DOUBLE ) {
-        for (int d=0; d<ndim; d++)
-            memcpy(&coords(0,d),mesh->coords[d],N_nodes*sizeof(double));
-    } else {
-        ERROR("Unsupported format");
+    int mesh_type = mesh->datatype;
+    for (int d=0; d<ndim; d++) {
+        Array<TYPE> data2( N_nodes );
+        copyData<TYPE>( data2, mesh_type, mesh->coords[d] );
+        memcpy( &coords(0,d), data2.data(), N_nodes*sizeof(TYPE) );
     }
     auto zones = mesh->zones;
     int N_zones = zones->nzones;
     int ndim_zones = zones->ndims;
     ASSERT( zones->nshapes==1 );
-    int type = zones->shapetype[0];
+    int shape_type = zones->shapetype[0];
     int shapesize = zones->shapesize[0];
     tri.resize(N_zones,shapesize);
     for (int i=0; i<N_zones; i++) {
@@ -335,22 +377,23 @@ void readTriMesh( DBfile* fid, const std::string& meshname, Array<double>& coord
     }
     DBFreeUcdmesh( mesh );
 }
+template<class TYPE>
 void writeTriMeshVariable( DBfile* fid, int ndim, const std::string& meshname,
-    const std::string& varname, const Array<double>& data, VariableType type )
+    const std::string& varname, const Array<TYPE>& data, VariableType type )
 {
     int nvars = 0;
     int vartype = 0;
-    const double *vars[10] = { nullptr };
-    if ( type == NodeVariable ) {
+    const TYPE *vars[10] = { nullptr };
+    if ( type == VariableType::NodeVariable ) {
         vartype = DB_NODECENT;
         nvars = data.size(1);
         for (int i=0; i<nvars; i++)
             vars[i] = &data(0,i);
-    } else if ( type == EdgeVariable ) {
+    } else if ( type == VariableType::EdgeVariable ) {
         ERROR("Not finished");
-    } else if ( type == SurfaceVariable ) {
+    } else if ( type == VariableType::SurfaceVariable ) {
         ERROR("Not finished");
-    } else if ( type == VolumeVariable ) {
+    } else if ( type == VariableType::VolumeVariable ) {
         vartype = DB_ZONECENT;
         nvars = data.size(1);
         for (int i=0; i<nvars; i++)
@@ -366,18 +409,19 @@ void writeTriMeshVariable( DBfile* fid, int ndim, const std::string& meshname,
     for (int i=0; i<nvars; i++)
         varnames[i] = const_cast<char*>(var_names[i].c_str());
     DBPutUcdvar( fid, varname.c_str(), meshname.c_str(), nvars,
-        varnames.data(), vars, data.size(0), nullptr, 0, DB_DOUBLE, vartype, nullptr );
+        varnames.data(), vars, data.size(0), nullptr, 0, getType<TYPE>(), vartype, nullptr );
 }
-Array<double> readTriMeshVariable( DBfile* fid, const std::string& varname )
+template<class TYPE>
+Array<TYPE> readTriMeshVariable( DBfile* fid, const std::string& varname )
 {
     auto var = DBGetUcdvar( fid, varname.c_str() );
     ASSERT( var != nullptr );
-    Array<double> data( var->nels, var->nvals );
-    if ( var->datatype == DB_DOUBLE ) {
-        for (int i=0; i<var->nvals; i++)
-            memcpy( &data(0,i), var->vals[i], var->nels*sizeof(double) );
-    } else {
-        ERROR("Unsupported format");
+    Array<TYPE> data( var->nels, var->nvals );
+    int type = var->datatype;
+    for (int i=0; i<var->nvals; i++) {
+        Array<TYPE> data2( var->nels );
+        copyData<TYPE>( data2, type, var->vals[i] );
+        memcpy( &data(0,i), data2.data(), var->nels*sizeof(TYPE) );
     }
     DBFreeUcdvar( var );
     return data;
@@ -420,18 +464,39 @@ void writeMultiVar( DBfile* fid, const std::string& varname,
 
 
 // Explicit instantiations
-template void silo::write<int>(    DBfile*, const std::string&, const std::vector<int>&    );
-template void silo::write<float>(  DBfile*, const std::string&, const std::vector<float>&  );
-template void silo::write<double>( DBfile*, const std::string&, const std::vector<double>& );
 template std::vector<int>    silo::read<int>(    DBfile* fid, const std::string& varname );
 template std::vector<float>  silo::read<float>(  DBfile* fid, const std::string& varname );
 template std::vector<double> silo::read<double>( DBfile* fid, const std::string& varname );
 template void silo::writeUniformMesh<1>( DBfile*, const std::string&, const std::array<double,2>&, const std::array<int,1>& );
 template void silo::writeUniformMesh<2>( DBfile*, const std::string&, const std::array<double,4>&, const std::array<int,2>& );
 template void silo::writeUniformMesh<3>( DBfile*, const std::string&, const std::array<double,6>&, const std::array<int,3>& );
-template void silo::writeUniformMeshVariable<1>( DBfile*, const std::string&, const std::array<int,1>&, const std::string&, const Array<double>&, silo::VariableType );
-template void silo::writeUniformMeshVariable<2>( DBfile*, const std::string&, const std::array<int,2>&, const std::string&, const Array<double>&, silo::VariableType );
-template void silo::writeUniformMeshVariable<3>( DBfile*, const std::string&, const std::array<int,3>&, const std::string&, const Array<double>&, silo::VariableType );
+
+namespace {
+    template<class TYPE>
+    static void siloClassInstantiator() {
+        // use **all** functions here so that they get instantiated
+        const auto type = silo::VariableType::NullVariable;
+        Array<TYPE> tmp;
+        Array<int> tri;
+        silo::write<TYPE>( nullptr, "", std::vector<TYPE>() );
+        silo::read<TYPE>( nullptr, "" );
+        silo::writeUniformMeshVariable<1,TYPE>( nullptr, "", std::array<int,1>(), "", Array<TYPE>(), type );
+        silo::writeUniformMeshVariable<2,TYPE>( nullptr, "", std::array<int,2>(), "", Array<TYPE>(), type );
+        silo::writeUniformMeshVariable<3,TYPE>( nullptr, "", std::array<int,3>(), "", Array<TYPE>(), type );
+        silo::writePointMesh( nullptr, "", 0, 0, (const TYPE**) nullptr );
+        silo::writePointMeshVariable<TYPE>( nullptr, "", "", Array<TYPE>() );
+        silo::writeTriMesh<TYPE>( nullptr, "", 0, 0, 0, (const TYPE**) nullptr, 0, (const int**) nullptr );
+        silo::writeTriMeshVariable<TYPE>( nullptr, 0, "", "", Array<TYPE>(), type );
+        silo::readUniformMeshVariable<TYPE>( nullptr, "" );
+        silo::readPointMesh<TYPE>( nullptr, "" );
+        silo::readPointMeshVariable<TYPE>( nullptr, "" );
+        silo::readTriMesh<TYPE>( nullptr, "", tmp, tri );
+        silo::readTriMeshVariable<TYPE>( nullptr, "" );
+    };
+    template void siloClassInstantiator<double>();
+    template void siloClassInstantiator<float>();
+    template void siloClassInstantiator<int>();
+}
 
 #else
 
diff --git a/IO/silo.h b/IO/silo.h
index c8e73b71..29e68245 100644
--- a/IO/silo.h
+++ b/IO/silo.h
@@ -23,7 +23,9 @@ namespace silo {
 
 enum FileMode { READ, WRITE, CREATE };
 
-enum VariableType { NodeVariable=1, EdgeVariable=2, SurfaceVariable=2, VolumeVariable=3, NullVariable=0 };
+enum class VariableType : int { NodeVariable=1, EdgeVariable=2, SurfaceVariable=2, VolumeVariable=3, NullVariable=0 };
+
+enum class VariableDataType { DOUBLE, FLOAT, INT, UNKNOWN };
 
 
 /*!
@@ -44,6 +46,15 @@ DBfile* open( const std::string& filename, FileMode mode );
 void close( DBfile* fid );
 
 
+/*!
+ * @brief  Get the variable type
+ * @detailed  This function returns the type of variable data
+ * @param[in] fid           Handle to the open file
+ * @param[in] name          Name of variable
+*/
+VariableDataType varDataType( DBfile *dbfile, const std::string& name );
+
+
 /*!
  * @brief  Write data to silo
  * @detailed  This function writes an arbitrary array to silo
@@ -101,9 +112,9 @@ void readUniformMesh( DBfile* fid, const std::string& meshname,
  * @param[in] data          Variable data
  * @param[in] type          Variable type
 */
-template<int NDIM>
+template< int NDIM, class TYPE >
 void writeUniformMeshVariable( DBfile* fid, const std::string& meshname, const std::array<int,NDIM>& N,
-    const std::string& varname, const Array<double>& data, VariableType type );
+    const std::string& varname, const Array<TYPE>& data, VariableType type );
 
 
 /*!
@@ -113,7 +124,8 @@ void writeUniformMeshVariable( DBfile* fid, const std::string& meshname, const s
  * @param[in] varname       Variable name
  * @return                  Variable data
 */
-Array<double> readUniformMeshVariable( DBfile* fid, const std::string& varname );
+template<class TYPE>
+Array<TYPE> readUniformMeshVariable( DBfile* fid, const std::string& varname );
 
 
 /*!
@@ -125,8 +137,9 @@ Array<double> readUniformMeshVariable( DBfile* fid, const std::string& varname )
  * @param[in] N             Number of points
  * @param[in] coords        Coordinates of the points
 */
+template<class TYPE>
 void writePointMesh( DBfile* fid, const std::string& meshname,
-    int ndim, int N, const double *coords[] );
+    int ndim, int N, const TYPE *coords[] );
 
 
 /*!
@@ -136,7 +149,8 @@ void writePointMesh( DBfile* fid, const std::string& meshname,
  * @param[in] meshname      Mesh name
  * @return                  Returns the coordinates as a N x ndim array 
 */
-Array<double> readPointMesh( DBfile* fid, const std::string& meshname );
+template<class TYPE>
+Array<TYPE> readPointMesh( DBfile* fid, const std::string& meshname );
 
 
 /*!
@@ -147,8 +161,9 @@ Array<double> readPointMesh( DBfile* fid, const std::string& meshname );
  * @param[in] varname       Variable name
  * @param[in] data          Variable data
 */
+template<class TYPE>
 void writePointMeshVariable( DBfile* fid, const std::string& meshname,
-    const std::string& varname, const Array<double>& data );
+    const std::string& varname, const Array<TYPE>& data );
 
 
 /*!
@@ -158,7 +173,8 @@ void writePointMeshVariable( DBfile* fid, const std::string& meshname,
  * @param[in] varname       Variable name
  * @return                  Variable data
 */
-Array<double> readPointMeshVariable( DBfile* fid, const std::string& varname );
+template<class TYPE>
+Array<TYPE> readPointMeshVariable( DBfile* fid, const std::string& varname );
 
 
 /*!
@@ -173,8 +189,9 @@ Array<double> readPointMeshVariable( DBfile* fid, const std::string& varname );
  * @param[in] N_tri         Number of triangles
  * @param[in] tri           Coordinates of the points
 */
+template<class TYPE>
 void writeTriMesh( DBfile* fid, const std::string& meshname,
-    int ndim, int ndim_tri, int N, const double *coords[], int N_tri, const int *tri[] );
+    int ndim, int ndim_tri, int N, const TYPE *coords[], int N_tri, const int *tri[] );
 
 
 /*!
@@ -185,7 +202,8 @@ void writeTriMesh( DBfile* fid, const std::string& meshname,
  * @param[in] coords        Coordinates of the points
  * @param[in] tri           Coordinates of the points
 */
-void readTriMesh( DBfile* fid, const std::string& meshname, Array<double>& coords, Array<int>& tri );
+template<class TYPE>
+void readTriMesh( DBfile* fid, const std::string& meshname, Array<TYPE>& coords, Array<int>& tri );
 
 
 /*!
@@ -195,9 +213,11 @@ void readTriMesh( DBfile* fid, const std::string& meshname, Array<double>& coord
  * @param[in] meshname      Mesh name
  * @param[in] varname       Variable name
  * @param[in] data          Variable data
+ * @param[in] type          Variable type
 */
+template<class TYPE>
 void writeTriMeshVariable( DBfile* fid, int ndim, const std::string& meshname,
-    const std::string& varname, const Array<double>& data, VariableType type );
+    const std::string& varname, const Array<TYPE>& data, VariableType type );
 
 
 /*!
@@ -207,7 +227,8 @@ void writeTriMeshVariable( DBfile* fid, int ndim, const std::string& meshname,
  * @param[in] varname       Variable name
  * @return                  Variable data
 */
-Array<double> readTriMeshVariable( DBfile* fid, const std::string& varname );
+template<class TYPE>
+Array<TYPE> readTriMeshVariable( DBfile* fid, const std::string& varname );
 
 
 /*!
diff --git a/tests/ComponentLabel.cpp b/tests/ComponentLabel.cpp
index df3c977a..067c40e4 100644
--- a/tests/ComponentLabel.cpp
+++ b/tests/ComponentLabel.cpp
@@ -442,31 +442,31 @@ int main(int argc, char **argv)
     std::shared_ptr<IO::Variable> PhaseIDVar( new IO::Variable() );
 
     PhaseVar->name = "phase";
-    PhaseVar->type = IO::VolumeVariable;
+    PhaseVar->type = IO::VariableType::VolumeVariable;
     PhaseVar->dim = 1;
     PhaseVar->data.resize(Nx-2,Ny-2,Nz-2);
     meshData[0].vars.push_back(PhaseVar);
 
     SignDistVar->name = "SignDist";
-    SignDistVar->type = IO::VolumeVariable;
+    SignDistVar->type = IO::VariableType::VolumeVariable;
     SignDistVar->dim = 1;
     SignDistVar->data.resize(Nx-2,Ny-2,Nz-2);
     meshData[0].vars.push_back(SignDistVar);
 
     LabelNWPVar->name = "LabelNWP";
-    LabelNWPVar->type = IO::VolumeVariable;
+    LabelNWPVar->type = IO::VariableType::VolumeVariable;
     LabelNWPVar->dim = 1;
     LabelNWPVar->data.resize(Nx-2,Ny-2,Nz-2);
     meshData[0].vars.push_back(LabelNWPVar);
 
     LabelWPVar->name = "LabelWP";
-    LabelWPVar->type = IO::VolumeVariable;
+    LabelWPVar->type = IO::VariableType::VolumeVariable;
     LabelWPVar->dim = 1;
     LabelWPVar->data.resize(Nx-2,Ny-2,Nz-2);
     meshData[0].vars.push_back(LabelWPVar);
 
     PhaseIDVar->name = "PhaseID";
-    PhaseIDVar->type = IO::VolumeVariable;
+    PhaseIDVar->type = IO::VariableType::VolumeVariable;
     PhaseIDVar->dim = 1;
     PhaseIDVar->data.resize(Nx-2,Ny-2,Nz-2);
     meshData[0].vars.push_back(PhaseIDVar);
diff --git a/tests/TestBlobIdentify.cpp b/tests/TestBlobIdentify.cpp
index bf588797..fbbae858 100644
--- a/tests/TestBlobIdentify.cpp
+++ b/tests/TestBlobIdentify.cpp
@@ -210,17 +210,17 @@ int main(int argc, char **argv)
     std::shared_ptr<IO::Variable> SignDistVar( new IO::Variable() );
     std::shared_ptr<IO::Variable> BlobIDVar( new IO::Variable() );
     PhaseVar->name = "phase";
-    PhaseVar->type = IO::VolumeVariable;
+    PhaseVar->type = IO::VariableType::VolumeVariable;
     PhaseVar->dim = 1;
     PhaseVar->data.resize(nx,ny,nz);
     meshData[0].vars.push_back(PhaseVar);
     SignDistVar->name = "SignDist";
-    SignDistVar->type = IO::VolumeVariable;
+    SignDistVar->type = IO::VariableType::VolumeVariable;
     SignDistVar->dim = 1;
     SignDistVar->data.resize(nx,ny,nz);
     meshData[0].vars.push_back(SignDistVar);
     BlobIDVar->name = "BlobID";
-    BlobIDVar->type = IO::VolumeVariable;
+    BlobIDVar->type = IO::VariableType::VolumeVariable;
     BlobIDVar->dim = 1;
     BlobIDVar->data.resize(nx,ny,nz);
     meshData[0].vars.push_back(BlobIDVar);
diff --git a/tests/TestBubble.cpp b/tests/TestBubble.cpp
index 840b4313..727e2c67 100644
--- a/tests/TestBubble.cpp
+++ b/tests/TestBubble.cpp
@@ -2175,7 +2175,7 @@ int main(int argc, char **argv)
                 shared_ptr<IO::Variable> dist( new IO::Variable() );
                 dist->name = "distance";
                 dist->dim = 1;
-                dist->type = IO::NodeVariable;
+                dist->type = IO::VariableType::NodeVariable;
                 dist->data.resize(3,mesh->A.size());
                 for (size_t i=0; i<mesh->A.size(); i++) {
                     const Point& a = mesh->A[i];
diff --git a/tests/TestWriter.cpp b/tests/TestWriter.cpp
index 52377226..f46dcaf0 100644
--- a/tests/TestWriter.cpp
+++ b/tests/TestWriter.cpp
@@ -16,12 +16,12 @@
 
 inline bool approx_equal( const Point& A, const Point& B )
 {
-    double tol = 1e-8*sqrt(A.x*A.x+A.y*A.y+A.z*A.z); 
+    double tol = 1e-7*sqrt(A.x*A.x+A.y*A.y+A.z*A.z); 
     return fabs(A.x-B.x)<=tol && fabs(A.y-B.y)<=tol && fabs(A.z-B.z)<=tol;
 }
 inline bool approx_equal( const double& A, const double& B )
 {
-    return fabs(A-B)<=1e-8*fabs(A+B);
+    return fabs(A-B)<=1e-7*fabs(A+B);
 }
 
 
@@ -32,7 +32,7 @@ inline double distance( const Point& p )
 
 
 // Test writing and reading the given format
-void testWriter( const std::string& format, const std::vector<IO::MeshDataStruct>& meshData, UnitTest& ut )
+void testWriter( const std::string& format, std::vector<IO::MeshDataStruct>& meshData, UnitTest& ut )
 {
     int rank, nprocs;
     MPI_Comm comm = MPI_COMM_WORLD;
@@ -40,9 +40,27 @@ void testWriter( const std::string& format, const std::vector<IO::MeshDataStruct
     MPI_Comm_size(comm,&nprocs);
     MPI_Barrier(comm);
 
+    // Get the format
+    std::string format2 = format;
+    auto precision = IO::DataType::Double;
+    if ( format == "silo-double" ) {
+        format2 = "silo";
+        precision = IO::DataType::Double;
+    } else if ( format == "silo-float" ) {
+        format2 = "silo";
+        precision = IO::DataType::Float;
+    }
+
+    // Set the precision for the variables
+    for ( auto& data : meshData ) {
+        data.precision = precision;
+        for ( auto& var : data.vars )
+            var->precision = precision;
+    }
+
     // Write the data
     PROFILE_START(format+"-write");
-    IO::initialize( "test_"+format, format, false );
+    IO::initialize( "test_"+format, format2, false );
     IO::writeData( 0, meshData, comm );
     IO::writeData( 3, meshData, comm );
     MPI_Barrier(comm);
@@ -53,7 +71,7 @@ void testWriter( const std::string& format, const std::vector<IO::MeshDataStruct
     std::string summary_name;
     if ( format=="old" || format=="new" )
         summary_name = "summary.LBM";
-    else if ( format == "silo" )
+    else if ( format=="silo-float" || format=="silo-double" )
         summary_name = "LBM.visit";
     else
         ERROR("Unknown format");
@@ -70,9 +88,9 @@ void testWriter( const std::string& format, const std::vector<IO::MeshDataStruct
     std::vector<std::string> timesteps = IO::readTimesteps( path + "/" + summary_name );
     PROFILE_STOP(format+"-read-timesteps");
     if ( timesteps.size()==2 )
-        ut.passes("Corrent number of timesteps");
+        ut.passes(format+": Corrent number of timesteps");
     else
-        ut.failure("Incorrent number of timesteps");
+        ut.failure(format+": Incorrent number of timesteps");
 
     // Check the mesh lists
     for ( const auto& timestep : timesteps ) {
@@ -81,17 +99,17 @@ void testWriter( const std::string& format, const std::vector<IO::MeshDataStruct
         auto databaseList = IO::getMeshList(path,timestep);
         PROFILE_STOP(format+"-read-getMeshList");
         if ( databaseList.size()==meshData.size() )
-            ut.passes("Corrent number of meshes found");
+            ut.passes(format+": Corrent number of meshes found");
         else
-            ut.failure("Incorrent number of meshes found");
+            ut.failure(format+": Incorrent number of meshes found");
         // Check the number of domains for each mesh
         bool pass = true;
         for ( const auto& database : databaseList )
             pass = pass && (int)database.domains.size()==nprocs;
         if ( pass ) {
-            ut.passes("Corrent number of domains for mesh");
+            ut.passes(format+": Corrent number of domains for mesh");
         } else {
-            ut.failure("Incorrent number of domains for mesh");
+            ut.failure(format+": Incorrent number of domains for mesh");
             continue;
         }
         // For each domain, load the mesh and check its data
@@ -161,9 +179,9 @@ void testWriter( const std::string& format, const std::vector<IO::MeshDataStruct
                 }
             }
             if ( pass ) {
-                ut.passes("Mesh \"" + database.name + "\" loaded correctly");
+                ut.passes(format+": Mesh \"" + database.name + "\" loaded correctly");
             } else {
-                ut.failure("Mesh \"" + database.name + "\" did not load correctly");
+                ut.failure(format+": Mesh \"" + database.name + "\" did not load correctly");
                 continue;
             }
             // Load the variables and check their data
@@ -196,9 +214,9 @@ void testWriter( const std::string& format, const std::vector<IO::MeshDataStruct
                             pass = pass && approx_equal(var1.data(m),var2.data(m));
                     }
                     if ( pass ) {
-                        ut.passes("Variable \"" + variable->name + "\" matched");
+                        ut.passes(format+": Variable \"" + variable->name + "\" matched");
                     } else {
-                        ut.failure("Variable \"" + variable->name + "\" did not match");
+                        ut.failure(format+": Variable \"" + variable->name + "\" did not match");
                         break;
                     }
                 }
@@ -260,20 +278,22 @@ int main(int argc, char **argv)
     std::shared_ptr<IO::DomainMesh> domain( new IO::DomainMesh(rank_data,6,7,8,1.0,1.0,1.0) );
 
     // Create the variables
-    std::shared_ptr<IO::Variable> set_node_mag( new IO::Variable(1,IO::NodeVariable,"Node_set_mag") );
-    std::shared_ptr<IO::Variable> set_node_vec( new IO::Variable(3,IO::NodeVariable,"Node_set_vec") );
-    std::shared_ptr<IO::Variable> list_node_mag( new IO::Variable(1,IO::NodeVariable,"Node_list_mag") );
-    std::shared_ptr<IO::Variable> list_node_vec( new IO::Variable(3,IO::NodeVariable,"Node_list_vec") );
-    std::shared_ptr<IO::Variable> point_node_mag( new IO::Variable(1,IO::NodeVariable,"Node_point_mag") );
-    std::shared_ptr<IO::Variable> point_node_vec( new IO::Variable(3,IO::NodeVariable,"Node_point_vec") );
-    std::shared_ptr<IO::Variable> domain_node_mag( new IO::Variable(1,IO::NodeVariable,"Node_domain_mag") );
-    std::shared_ptr<IO::Variable> domain_node_vec( new IO::Variable(3,IO::NodeVariable,"Node_domain_vec") );
-    std::shared_ptr<IO::Variable> set_cell_mag( new IO::Variable(1,IO::VolumeVariable,"Cell_set_mag") );
-    std::shared_ptr<IO::Variable> set_cell_vec( new IO::Variable(3,IO::VolumeVariable,"Cell_set_vec") );
-    std::shared_ptr<IO::Variable> list_cell_mag( new IO::Variable(1,IO::VolumeVariable,"Cell_list_mag") );
-    std::shared_ptr<IO::Variable> list_cell_vec( new IO::Variable(3,IO::VolumeVariable,"Cell_list_vec") );
-    std::shared_ptr<IO::Variable> domain_cell_mag( new IO::Variable(1,IO::VolumeVariable,"Cell_domain_mag") );
-    std::shared_ptr<IO::Variable> domain_cell_vec( new IO::Variable(3,IO::VolumeVariable,"Cell_domain_vec") );
+    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") );
     point_node_mag->data.resize( N_points );
     point_node_vec->data.resize( N_points, 3 );
     for (int i=0; i<N_points; i++) {
@@ -363,7 +383,8 @@ int main(int argc, char **argv)
     // Run the tests
     testWriter( "old", meshData, ut );
     testWriter( "new", meshData, ut );
-    testWriter( "silo", meshData, ut );
+    testWriter( "silo-double", meshData, ut );
+    testWriter( "silo-float", meshData, ut );
 
     // Finished
     ut.report();
diff --git a/tests/lbpm_color_macro_simulator.cpp b/tests/lbpm_color_macro_simulator.cpp
index 2422bda7..3d718fb1 100644
--- a/tests/lbpm_color_macro_simulator.cpp
+++ b/tests/lbpm_color_macro_simulator.cpp
@@ -714,22 +714,22 @@ int main(int argc, char **argv)
     std::shared_ptr<IO::Variable> SignDistVar( new IO::Variable() );
     std::shared_ptr<IO::Variable> BlobIDVar( new IO::Variable() );
     PhaseVar->name = "phase";
-    PhaseVar->type = IO::VolumeVariable;
+    PhaseVar->type = IO::VariableType::VolumeVariable;
     PhaseVar->dim = 1;
     PhaseVar->data.resize(Nx-2,Ny-2,Nz-2);
     meshData[0].vars.push_back(PhaseVar);
     PressVar->name = "Pressure";
-    PressVar->type = IO::VolumeVariable;
+    PressVar->type = IO::VariableType::VolumeVariable;
     PressVar->dim = 1;
     PressVar->data.resize(Nx-2,Ny-2,Nz-2);
     meshData[0].vars.push_back(PressVar);
     SignDistVar->name = "SignDist";
-    SignDistVar->type = IO::VolumeVariable;
+    SignDistVar->type = IO::VariableType::VolumeVariable;
     SignDistVar->dim = 1;
     SignDistVar->data.resize(Nx-2,Ny-2,Nz-2);
     meshData[0].vars.push_back(SignDistVar);
     BlobIDVar->name = "BlobID";
-    BlobIDVar->type = IO::VolumeVariable;
+    BlobIDVar->type = IO::VariableType::VolumeVariable;
     BlobIDVar->dim = 1;
     BlobIDVar->data.resize(Nx-2,Ny-2,Nz-2);
     meshData[0].vars.push_back(BlobIDVar);
diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp
index d799de62..639b1cd5 100644
--- a/tests/lbpm_color_simulator.cpp
+++ b/tests/lbpm_color_simulator.cpp
@@ -755,22 +755,22 @@ int main(int argc, char **argv)
     std::shared_ptr<IO::Variable> SignDistVar( new IO::Variable() );
     std::shared_ptr<IO::Variable> BlobIDVar( new IO::Variable() );
     PhaseVar->name = "phase";
-    PhaseVar->type = IO::VolumeVariable;
+    PhaseVar->type = IO::VariableType::VolumeVariable;
     PhaseVar->dim = 1;
     PhaseVar->data.resize(Nx-2,Ny-2,Nz-2);
     meshData[0].vars.push_back(PhaseVar);
     PressVar->name = "Pressure";
-    PressVar->type = IO::VolumeVariable;
+    PressVar->type = IO::VariableType::VolumeVariable;
     PressVar->dim = 1;
     PressVar->data.resize(Nx-2,Ny-2,Nz-2);
     meshData[0].vars.push_back(PressVar);
     SignDistVar->name = "SignDist";
-    SignDistVar->type = IO::VolumeVariable;
+    SignDistVar->type = IO::VariableType::VolumeVariable;
     SignDistVar->dim = 1;
     SignDistVar->data.resize(Nx-2,Ny-2,Nz-2);
     meshData[0].vars.push_back(SignDistVar);
     BlobIDVar->name = "BlobID";
-    BlobIDVar->type = IO::VolumeVariable;
+    BlobIDVar->type = IO::VariableType::VolumeVariable;
     BlobIDVar->dim = 1;
     BlobIDVar->data.resize(Nx-2,Ny-2,Nz-2);
     meshData[0].vars.push_back(BlobIDVar);
diff --git a/tests/lbpm_uCT_maskfilter.cpp b/tests/lbpm_uCT_maskfilter.cpp
index f5ba65d9..22bcb45f 100644
--- a/tests/lbpm_uCT_maskfilter.cpp
+++ b/tests/lbpm_uCT_maskfilter.cpp
@@ -329,7 +329,7 @@ int main(int argc, char **argv)
         // Source data
         std::shared_ptr<IO::Variable> OrigData( new IO::Variable() );
 	    OrigData->name = "Source Data";
-	    OrigData->type = IO::VolumeVariable;
+	    OrigData->type = IO::VariableType::VolumeVariable;
 	    OrigData->dim = 1;
 	    OrigData->data.resize(Nx[i],Ny[i],Nz[i]);
 	    meshData[i].vars.push_back(OrigData);
@@ -337,7 +337,7 @@ int main(int argc, char **argv)
         // Non-Local Mean
         std::shared_ptr<IO::Variable> NonLocMean( new IO::Variable() );
 	    NonLocMean->name = "Non-Local Mean";
-	    NonLocMean->type = IO::VolumeVariable;
+	    NonLocMean->type = IO::VariableType::VolumeVariable;
 	    NonLocMean->dim = 1;
 	    NonLocMean->data.resize(Nx[i],Ny[i],Nz[i]);
 	    meshData[i].vars.push_back(NonLocMean);
@@ -345,7 +345,7 @@ int main(int argc, char **argv)
         // Segmented Data
         std::shared_ptr<IO::Variable> SegData( new IO::Variable() );
 	    SegData->name = "Segmented Data";
-	    SegData->type = IO::VolumeVariable;
+	    SegData->type = IO::VariableType::VolumeVariable;
 	    SegData->dim = 1;
 	    SegData->data.resize(Nx[i],Ny[i],Nz[i]);
 	    meshData[i].vars.push_back(SegData);
@@ -353,7 +353,7 @@ int main(int argc, char **argv)
         // Signed Distance
         std::shared_ptr<IO::Variable> DistData( new IO::Variable() );
 	    DistData->name = "Signed Distance";
-	    DistData->type = IO::VolumeVariable;
+	    DistData->type = IO::VariableType::VolumeVariable;
 	    DistData->dim = 1;
 	    DistData->data.resize(Nx[i],Ny[i],Nz[i]);
 	    meshData[i].vars.push_back(DistData);
@@ -361,7 +361,7 @@ int main(int argc, char **argv)
         // Smoothed Data
         std::shared_ptr<IO::Variable> SmoothData( new IO::Variable() );
 	    SmoothData->name = "Smoothed Data";
-	    SmoothData->type = IO::VolumeVariable;
+	    SmoothData->type = IO::VariableType::VolumeVariable;
 	    SmoothData->dim = 1;
 	    SmoothData->data.resize(Nx[i],Ny[i],Nz[i]);
 	    meshData[i].vars.push_back(SmoothData);
@@ -370,21 +370,21 @@ int main(int argc, char **argv)
     #if 0
         std::shared_ptr<IO::Variable> filter_Mean_var( new IO::Variable() );
 	    filter_Mean_var->name = "Mean";
-	    filter_Mean_var->type = IO::VolumeVariable;
+	    filter_Mean_var->type = IO::VariableType::VolumeVariable;
 	    filter_Mean_var->dim = 1;
 	    filter_Mean_var->data.resize(Nx[0],Ny[0],Nz[0]);
 	    meshData[0].vars.push_back(filter_Mean_var);
         fillDouble[0]->copy( filter_Mean, filter_Mean_var->data );
         std::shared_ptr<IO::Variable> filter_Dist1_var( new IO::Variable() );
 	    filter_Dist1_var->name = "Dist1";
-	    filter_Dist1_var->type = IO::VolumeVariable;
+	    filter_Dist1_var->type = IO::VariableType::VolumeVariable;
 	    filter_Dist1_var->dim = 1;
 	    filter_Dist1_var->data.resize(Nx[0],Ny[0],Nz[0]);
 	    meshData[0].vars.push_back(filter_Dist1_var);
         fillDouble[0]->copy( filter_Dist1, filter_Dist1_var->data );
         std::shared_ptr<IO::Variable> filter_Dist2_var( new IO::Variable() );
 	    filter_Dist2_var->name = "Dist2";
-	    filter_Dist2_var->type = IO::VolumeVariable;
+	    filter_Dist2_var->type = IO::VariableType::VolumeVariable;
 	    filter_Dist2_var->dim = 1;
 	    filter_Dist2_var->data.resize(Nx[0],Ny[0],Nz[0]);
 	    meshData[0].vars.push_back(filter_Dist2_var);
diff --git a/tests/lbpm_uCT_pp.cpp b/tests/lbpm_uCT_pp.cpp
index 85300436..6137a8b5 100644
--- a/tests/lbpm_uCT_pp.cpp
+++ b/tests/lbpm_uCT_pp.cpp
@@ -291,7 +291,7 @@ int main(int argc, char **argv)
         // Source data
         std::shared_ptr<IO::Variable> OrigData( new IO::Variable() );
 	    OrigData->name = "Source Data";
-	    OrigData->type = IO::VolumeVariable;
+	    OrigData->type = IO::VariableType::VolumeVariable;
 	    OrigData->dim = 1;
 	    OrigData->data.resize(Nx[i],Ny[i],Nz[i]);
 	    meshData[i].vars.push_back(OrigData);
@@ -299,7 +299,7 @@ int main(int argc, char **argv)
         // Non-Local Mean
         std::shared_ptr<IO::Variable> NonLocMean( new IO::Variable() );
 	    NonLocMean->name = "Non-Local Mean";
-	    NonLocMean->type = IO::VolumeVariable;
+	    NonLocMean->type = IO::VariableType::VolumeVariable;
 	    NonLocMean->dim = 1;
 	    NonLocMean->data.resize(Nx[i],Ny[i],Nz[i]);
 	    meshData[i].vars.push_back(NonLocMean);
@@ -307,7 +307,7 @@ int main(int argc, char **argv)
         // Segmented Data
         std::shared_ptr<IO::Variable> SegData( new IO::Variable() );
 	    SegData->name = "Segmented Data";
-	    SegData->type = IO::VolumeVariable;
+	    SegData->type = IO::VariableType::VolumeVariable;
 	    SegData->dim = 1;
 	    SegData->data.resize(Nx[i],Ny[i],Nz[i]);
 	    meshData[i].vars.push_back(SegData);
@@ -315,7 +315,7 @@ int main(int argc, char **argv)
         // Signed Distance
         std::shared_ptr<IO::Variable> DistData( new IO::Variable() );
 	    DistData->name = "Signed Distance";
-	    DistData->type = IO::VolumeVariable;
+	    DistData->type = IO::VariableType::VolumeVariable;
 	    DistData->dim = 1;
 	    DistData->data.resize(Nx[i],Ny[i],Nz[i]);
 	    meshData[i].vars.push_back(DistData);
@@ -323,7 +323,7 @@ int main(int argc, char **argv)
         // Smoothed Data
         std::shared_ptr<IO::Variable> SmoothData( new IO::Variable() );
 	    SmoothData->name = "Smoothed Data";
-	    SmoothData->type = IO::VolumeVariable;
+	    SmoothData->type = IO::VariableType::VolumeVariable;
 	    SmoothData->dim = 1;
 	    SmoothData->data.resize(Nx[i],Ny[i],Nz[i]);
 	    meshData[i].vars.push_back(SmoothData);
@@ -332,21 +332,21 @@ int main(int argc, char **argv)
     #if 0
         std::shared_ptr<IO::Variable> filter_Mean_var( new IO::Variable() );
 	    filter_Mean_var->name = "Mean";
-	    filter_Mean_var->type = IO::VolumeVariable;
+	    filter_Mean_var->type = IO::VariableType::VolumeVariable;
 	    filter_Mean_var->dim = 1;
 	    filter_Mean_var->data.resize(Nx[0],Ny[0],Nz[0]);
 	    meshData[0].vars.push_back(filter_Mean_var);
         fillDouble[0]->copy( filter_Mean, filter_Mean_var->data );
         std::shared_ptr<IO::Variable> filter_Dist1_var( new IO::Variable() );
 	    filter_Dist1_var->name = "Dist1";
-	    filter_Dist1_var->type = IO::VolumeVariable;
+	    filter_Dist1_var->type = IO::VariableType::VolumeVariable;
 	    filter_Dist1_var->dim = 1;
 	    filter_Dist1_var->data.resize(Nx[0],Ny[0],Nz[0]);
 	    meshData[0].vars.push_back(filter_Dist1_var);
         fillDouble[0]->copy( filter_Dist1, filter_Dist1_var->data );
         std::shared_ptr<IO::Variable> filter_Dist2_var( new IO::Variable() );
 	    filter_Dist2_var->name = "Dist2";
-	    filter_Dist2_var->type = IO::VolumeVariable;
+	    filter_Dist2_var->type = IO::VariableType::VolumeVariable;
 	    filter_Dist2_var->dim = 1;
 	    filter_Dist2_var->data.resize(Nx[0],Ny[0],Nz[0]);
 	    meshData[0].vars.push_back(filter_Dist2_var);