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

Adding initial Silo writer. Modifying writer interface.

Browse Source
This commit is contained in:
Mark Berrill 2017-01-23 10:42:21 -05:00
parent 5329c1f659
commit e29faadcee
23 changed files with 983 additions and 229 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists.txt
View File

@ -106,6 +106,7 @@ IF ( NOT ONLY_BUILD_DOCS )
CONFIGURE_CUDA()
CONFIGURE_MIC()
CONFIGURE_NETCDF()
CONFIGURE_SILO()
CONFIGURE_LBPM()
CONFIGURE_TIMER( 0 "${${PROJ}_INSTALL_DIR}/null_timer" )
CONFIGURE_LINE_COVERAGE()

38
IO/MeshDatabase.cpp
View File

@ -238,42 +238,49 @@ void DatabaseEntry::read( const std::string& line )
// Gather the mesh databases from all processors
inline int tod( int N ) { return (N+7)/sizeof(double); }
std::vector<MeshDatabase> gatherAll( const std::vector<MeshDatabase>& meshes, MPI_Comm comm )
{
#ifdef USE_MPI
PROFILE_START("gatherAll");
PROFILE_START("gatherAll-pack",2);
int size = MPI_WORLD_SIZE();
// First pack the mesh data to local buffers
size_t localsize = 0;
int localsize = 0;
for (size_t i=0; i<meshes.size(); i++)
localsize += packsize(meshes[i]);
char *localbuf = new char[localsize];
size_t pos = 0;
localsize += tod(packsize(meshes[i]));
auto localbuf = new double[localsize];
int pos = 0;
for (size_t i=0; i<meshes.size(); i++) {
pack( meshes[i], &localbuf[pos] );
pos += packsize(meshes[i]);
pack( meshes[i], (char*) &localbuf[pos] );
pos += tod(packsize(meshes[i]));
}
PROFILE_STOP("gatherAll-pack",2);
// Get the number of bytes each processor will be sending/recieving
int sendsize = static_cast<int>(localsize);
int *recvsize = new int[size];
MPI_Allgather(&sendsize,1,MPI_INT,recvsize,1,MPI_INT,comm);
size_t globalsize = recvsize[0];
int *disp = new int[size];
PROFILE_START("gatherAll-send1",2);
auto recvsize = new int[size];
MPI_Allgather(&localsize,1,MPI_INT,recvsize,1,MPI_INT,comm);
int globalsize = recvsize[0];
auto disp = new int[size];
disp[0] = 0;
for (int i=1; i<size; i++) {
disp[i] = disp[i-1] + recvsize[i];
globalsize += recvsize[i];
}
PROFILE_STOP("gatherAll-send1",2);
// Send/recv the global data
char *globalbuf = new char[globalsize];
MPI_Allgatherv(localbuf,sendsize,MPI_CHAR,globalbuf,recvsize,disp,MPI_CHAR,comm);
PROFILE_START("gatherAll-send2",2);
auto globalbuf = new double[globalsize];
MPI_Allgatherv(localbuf,localsize,MPI_DOUBLE,globalbuf,recvsize,disp,MPI_DOUBLE,comm);
PROFILE_STOP("gatherAll-send2",2);
// Unpack the data
PROFILE_START("gatherAll-unpack",2);
std::map<std::string,MeshDatabase> data;
pos = 0;
while ( pos < globalsize ) {
MeshDatabase tmp;
unpack(tmp,&globalbuf[pos]);
pos += packsize(tmp);
unpack(tmp,(char*)&globalbuf[pos]);
pos += tod(packsize(tmp));
std::map<std::string,MeshDatabase>::iterator it = data.find(tmp.name);
if ( it==data.end() ) {
data[tmp.name] = tmp;
@ -300,6 +307,7 @@ std::vector<MeshDatabase> gatherAll( const std::vector<MeshDatabase>& meshes, MP
size_t i=0;
for (std::map<std::string,MeshDatabase>::iterator it=data.begin(); it!=data.end(); ++it, ++i)
data2[i] = it->second;
PROFILE_STOP("gatherAll-unpack",2);
PROFILE_STOP("gatherAll");
return data2;
#else

2
IO/MeshDatabase.h
View File

@ -56,7 +56,7 @@ struct MeshDatabase {
std::string name; //!< Name of the mesh
MeshType type; //!< Mesh type
std::string meshClass; //!< Mesh class
unsigned char format; //!< Data format
unsigned char format; //!< Data format (1: old, 2: new, 3: new (single), 4: silo)
std::vector<DatabaseEntry> domains; //!< List of the domains
std::vector<VariableDatabase> variables; //!< List of the variables
std::map<variable_id,DatabaseEntry> variable_data; //!< Data for the variables

53
IO/Reader.cpp
View File

@ -20,6 +20,18 @@ static inline void fgetl( char * str, int num, FILE * stream )
}
// Get the path to a file
std::string IO::getPath( const std::string& filename )
{
std::string file(filename);
size_t k1 = file.rfind(47);
size_t k2 = file.rfind(92);
if ( k1==std::string::npos ) { k1=0; }
if ( k2==std::string::npos ) { k2=0; }
return file.substr(0,std::max(k1,k2));
}
// List the timesteps in the given directors (dumps.LBPM)
std::vector<std::string> IO::readTimesteps( const std::string& filename )
{
@ -166,23 +178,50 @@ std::shared_ptr<IO::Variable> IO::getVariable( const std::string& path, const st
size_t N = atol(values[2].c_str());
size_t bytes = atol(values[3].c_str());
std::string precision = values[4];
char *data = new char[bytes];
size_t count = fread(data,1,bytes,fid);
fclose(fid);
ASSERT(count==bytes);
std::shared_ptr<IO::Variable> var( new IO::Variable() );
var->dim = dim;
var->type = static_cast<IO::VariableType>(type);
var->name = variable;
var->data.resize(N);
var->data.resize(N*dim);
if ( precision=="double" ) {
memcpy(var->data.data(),data,bytes);
size_t count = fread(var->data.data(),sizeof(double),N*dim,fid);
ASSERT(count*sizeof(double)==bytes);
} else {
ERROR("Format not implimented");
}
delete [] data;
fclose(fid);
return var;
}
/****************************************************
* Reformat the variable to match the mesh *
****************************************************/
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 ) {
ERROR("Not finished");
} else if ( var.type == EdgeVariable ) {
ERROR("Not finished");
} else if ( var.type == SurfaceVariable ) {
ERROR("Not finished");
} else if ( var.type == 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 } );
} else {
ERROR("Invalid variable type");
}
} else if ( mesh.className() == "PointList" ) {
ERROR("Not finished");
} else if ( mesh.className() == "TriMesh" ) {
ERROR("Not finished");
} else {
ERROR("Unknown mesh type");
}
}

26
IO/Reader.h
View File

@ -13,11 +13,15 @@
namespace IO {
//! Get the path to a file
std::string getPath( const std::string& filename );
//! List the timesteps in the given directors (dumps.LBPM)
std::vector<std::string> readTimesteps( const std::string& filename );
//! Read the list of variables for the given timestep
//! Read the list of mesh databases for the given timestep
std::vector<IO::MeshDatabase> getMeshList( const std::string& path, const std::string& timestep );
@ -26,11 +30,29 @@ std::shared_ptr<IO::Mesh> getMesh( const std::string& path, const std::string& t
const MeshDatabase& meshDatabase, int domain );
//! Read the given mesh domain
/*!
* @brief Read the given variable
* @details This function reads the variable data provided for the current timestep
* @param[in] path The path to the file
* @param[in] timestep The timestep iteration
* @param[in] meshDatabase The mesh database (see getMeshList)
* @param[in] domain The index of the domain we want to read
* @param[in] variable The variable name to read
* @return Returns the variable data as a linear array
*/
std::shared_ptr<IO::Variable> getVariable( const std::string& path, const std::string& timestep,
const MeshDatabase& meshDatabase, int domain, const std::string& variable );
/*!
* @brief Reformat the variable to match the mesh
* @details This function modifies the dimensions of the array to match the mesh
* @param[in] mesh The underlying mesh
* @param[in/out] variable The variable name to read
*/
void reformatVariable( const IO::Mesh& mesh, IO::Variable& var );
} // IO namespace
#endif

268
IO/Writer.cpp
View File

@ -1,6 +1,7 @@
#include "IO/Writer.h"
#include "IO/MeshDatabase.h"
#include "IO/IOHelpers.h"
#include "IO/silo.h"
#include "common/MPI_Helpers.h"
#include "common/Utilities.h"
#include "shared_ptr.h"
@ -11,11 +12,48 @@
#include <set>
static bool global_summary_created = false;
enum class Format { OLD, NEW, SILO, UNKNOWN };
/****************************************************
* Initialize the writer *
****************************************************/
static std::string global_IO_path;
static Format global_IO_format = Format::UNKNOWN;
void IO::initialize( const std::string& path, const std::string& format, bool append )
{
if ( path.empty() )
global_IO_path = ".";
else
global_IO_path = path;
if ( format == "old" )
global_IO_format = Format::OLD;
else if ( format == "new" )
global_IO_format = Format::NEW;
else if ( format == "silo" )
global_IO_format = Format::SILO;
else
ERROR("Unknown format");
int rank = comm_rank(MPI_COMM_WORLD);
if ( !append && rank==0 ) {
mkdir(path.c_str(),S_IRWXU|S_IRGRP);
std::string filename;
if ( global_IO_format==Format::OLD || global_IO_format==Format::NEW )
filename = global_IO_path + "/summary.LBM";
else if ( global_IO_format==Format::SILO )
filename = global_IO_path + "/LBM.visit";
else
ERROR("Unknown format");
auto fid = fopen(filename.c_str(),"wb");
fclose(fid);
}
}
// Write the mesh data in the original format
static std::vector<IO::MeshDatabase> writeMeshesOrigFormat( const std::vector<IO::MeshDataStruct>& meshData, const char* path )
static std::vector<IO::MeshDatabase> writeMeshesOrigFormat( const std::vector<IO::MeshDataStruct>& meshData, const std::string& path )
{
int rank = MPI_WORLD_RANK();
std::vector<IO::MeshDatabase> meshes_written;
@ -23,7 +61,7 @@ static std::vector<IO::MeshDatabase> writeMeshesOrigFormat( const std::vector<IO
char domainname[100], filename[100], fullpath[200];
sprintf(domainname,"%05i",rank);
sprintf(filename,"%s.%05i",meshData[i].meshName.c_str(),rank);
sprintf(fullpath,"%s/%s",path,filename);
sprintf(fullpath,"%s/%s",path.c_str(),filename);
FILE *fid = fopen(fullpath,"wb");
INSIST(fid!=NULL,std::string("Error opening file: ")+fullpath);
std::shared_ptr<IO::Mesh> mesh = meshData[i].mesh;
@ -79,14 +117,12 @@ static std::vector<IO::MeshDatabase> writeMeshesOrigFormat( const std::vector<IO
}
// Write a mesh (and variables) to a file
static IO::MeshDatabase write_domain( FILE *fid, const std::string& filename,
const IO::MeshDataStruct& mesh, int format )
// Create the database entry for the mesh data
static IO::MeshDatabase getDatabase( const std::string& filename, const IO::MeshDataStruct& mesh, int format )
{
int rank = MPI_WORLD_RANK();
char domainname[10];
sprintf(domainname,"%05i",rank);
int level = 0;
sprintf(domainname,"%s_%05i",mesh.meshName.c_str(),rank);
// Create the MeshDatabase
IO::MeshDatabase database;
database.name = mesh.meshName;
@ -97,8 +133,40 @@ static IO::MeshDatabase write_domain( FILE *fid, const std::string& filename,
IO::DatabaseEntry domain;
domain.name = domainname;
domain.file = filename;
domain.offset = ftell(fid);
domain.offset = -1;
database.domains.push_back(domain);
// Write the variables
for (size_t i=0; i<mesh.vars.size(); i++) {
// Add basic variable info
IO::VariableDatabase info;
info.name = mesh.vars[i]->name;
info.type = mesh.vars[i]->type;
info.dim = mesh.vars[i]->dim;
database.variables.push_back(info);
// Add domain variable info
IO::DatabaseEntry variable;
variable.name = mesh.vars[i]->name;
variable.file = filename;
variable.offset = -1;
std::pair<std::string,std::string> key(domain.name,mesh.vars[i]->name);
database.variable_data.insert(
std::pair<std::pair<std::string,std::string>,IO::DatabaseEntry>(key,variable) );
}
return database;
}
// Write a mesh (and variables) to a file
static IO::MeshDatabase write_domain( FILE *fid, const std::string& filename,
const IO::MeshDataStruct& mesh, int format )
{
const int level = 0;
int rank = MPI_WORLD_RANK();
// Create the MeshDatabase
IO::MeshDatabase database = getDatabase( filename, mesh, format );
// Write the mesh
IO::DatabaseEntry& domain = database.domains[0];
domain.offset = ftell(fid);
std::pair<size_t,void*> data = mesh.mesh->pack(level);
fprintf(fid,"Mesh: %s-%05i: %lu\n",mesh.meshName.c_str(),rank,data.first);
fwrite(data.second,1,data.first,fid);
@ -106,22 +174,12 @@ static IO::MeshDatabase write_domain( FILE *fid, const std::string& filename,
delete [] (char*) data.second;
// Write the variables
for (size_t i=0; i<mesh.vars.size(); i++) {
IO::DatabaseEntry variable;
variable.name = mesh.vars[i]->name;
variable.file = filename;
std::pair<std::string,std::string> key(domain.name,mesh.vars[i]->name);
IO::DatabaseEntry& variable = database.variable_data[key];
variable.offset = ftell(fid);
IO::VariableDatabase info;
info.name = variable.name;
info.type = mesh.vars[i]->type;
info.dim = mesh.vars[i]->dim;
database.variables.push_back(info);
std::pair<std::string,std::string> key(domainname,variable.name);
database.variable_data.insert(
std::pair<std::pair<std::string,std::string>,IO::DatabaseEntry>(key,variable) );
int dim = mesh.vars[i]->dim;
int type = static_cast<int>(mesh.vars[i]->type);
size_t N = mesh.vars[i]->data.length();
const double* data = N==0 ? NULL:mesh.vars[i]->data.data();
if ( type == static_cast<int>(IO::NullVariable) ) {
ERROR("Variable type not set");
}
@ -129,23 +187,101 @@ static IO::MeshDatabase write_domain( FILE *fid, const std::string& filename,
ASSERT(N==dim*N_mesh);
fprintf(fid,"Var: %s-%05i-%s: %i, %i, %lu, %lu, double\n",
database.name.c_str(), rank, variable.name.c_str(),
dim, type, N, dim*N*sizeof(double) );
fwrite(data,sizeof(double),dim*N,fid);
dim, type, N_mesh, N*sizeof(double) );
fwrite(mesh.vars[i]->data.data(),sizeof(double),N,fid);
fprintf(fid,"\n");
}
return database;
}
#ifdef USE_SILO
// Write a PointList mesh (and variables) to a file
static void writeSiloPointList( DBfile *fid, const IO::MeshDataStruct& meshData, IO::MeshDatabase database )
{
ERROR("Not finished yet");
}
// Write a TriMesh mesh (and variables) to a file
static void writeSiloTriMesh( DBfile *fid, const IO::MeshDataStruct& meshData, IO::MeshDatabase database )
{
ERROR("Not finished yet");
}
// Write a TriList mesh (and variables) to a file
static void writeSiloTriList( DBfile *fid, const IO::MeshDataStruct& meshData, IO::MeshDatabase database )
{
ERROR("Not finished yet");
}
// Write a DomainMesh mesh (and variables) to a file
static void writeSiloDomainMesh( DBfile *fid, const IO::MeshDataStruct& meshData, IO::MeshDatabase database )
{
const IO::DomainMesh& mesh = dynamic_cast<IO::DomainMesh&>( *meshData.mesh );
RankInfoStruct info( mesh.rank, mesh.nprocx, mesh.nprocy, mesh.nprocz );
std::array<double,6> range = { info.ix*mesh.Lx/info.nx, (info.ix+1)*mesh.Lx/info.nx,
info.jy*mesh.Ly/info.ny, (info.jy+1)*mesh.Ly/info.ny,
info.kz*mesh.Lz/info.nz, (info.kz+1)*mesh.Lz/info.nz };
std::array<int,3> N = { mesh.nx, mesh.ny, mesh.nz };
std::array<int,3> baseindex = { info.ix, info.jy, info.kz };
const std::string meshname = database.domains[0].name;
silo::writeUniformMesh<3>( fid, meshname, range, N );
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 );
}
}
// Write a mesh (and variables) to a file
static IO::MeshDatabase write_domain_silo( DBfile *fid, const std::string& filename,
const IO::MeshDataStruct& mesh, int format )
{
const int level = 0;
int rank = MPI_WORLD_RANK();
// Create the MeshDatabase
IO::MeshDatabase database = getDatabase( filename, mesh, format );
if ( database.meshClass=="PointList" ) {
writeSiloPointList( fid, mesh, database );
} else if ( database.meshClass=="TriMesh" ) {
writeSiloTriMesh( fid, mesh, database );
} else if ( database.meshClass=="TriList" ) {
writeSiloTriList( fid, mesh, database );
} else if ( database.meshClass=="DomainMesh" ) {
writeSiloDomainMesh( fid, mesh, database );
} else {
ERROR("Unknown mesh class");
}
return database;
}
// Write the summary file for silo
void writeSiloSummary( const std::vector<IO::MeshDatabase>& meshes_written, const std::string& filename )
{
auto fid = silo::open( filename, silo::CREATE );
for ( const auto& data : meshes_written ) {
std::vector<std::string> meshnames;
for ( const auto& tmp : data.domains )
meshnames.push_back( tmp.file + ":" + tmp.name );
std::vector<int> meshTypes(meshnames.size(),DB_QUAD_RECT); // Not correct for all types
silo::writeMultiMesh( fid, data.name, meshnames, meshTypes );
for (const auto& variable : data.variables ) {
std::vector<int> varTypes(meshnames.size(),DB_QUADVAR); // Not correct for all types
std::vector<std::string> varnames;
for ( const auto& tmp : data.domains )
varnames.push_back( tmp.file + ":" + variable.name );
silo::writeMultiVar( fid, variable.name, varnames, varTypes, 3, variable.dim );
}
}
silo::close( fid );
}
#endif
// Write the mesh data in the new format
static std::vector<IO::MeshDatabase> writeMeshesNewFormat(
const std::vector<IO::MeshDataStruct>& meshData, const char* path, int format )
const std::vector<IO::MeshDataStruct>& meshData, const std::string& path, int format )
{
int rank = MPI_WORLD_RANK();
std::vector<IO::MeshDatabase> meshes_written;
char filename[100], fullpath[200];
sprintf(filename,"%05i",rank);
sprintf(fullpath,"%s/%s",path,filename);
sprintf(fullpath,"%s/%s",path.c_str(),filename);
FILE *fid = fopen(fullpath,"wb");
for (size_t i=0; i<meshData.size(); i++) {
std::shared_ptr<IO::Mesh> mesh = meshData[i].mesh;
@ -156,25 +292,56 @@ static std::vector<IO::MeshDatabase> writeMeshesNewFormat(
}
// Write the mesh data
void IO::writeData( int timestep, const std::vector<IO::MeshDataStruct>& meshData, int format, MPI_Comm comm )
// Write the mesh data to silo
static std::vector<IO::MeshDatabase> writeMeshesSilo(
const std::vector<IO::MeshDataStruct>& meshData, const std::string& path, int format )
{
#ifdef USE_SILO
int rank = MPI_WORLD_RANK();
std::vector<IO::MeshDatabase> meshes_written;
char filename[100], fullpath[200];
sprintf(filename,"%05i.silo",rank);
sprintf(fullpath,"%s/%s",path.c_str(),filename);
auto fid = silo::open( fullpath, silo::CREATE );
for (size_t i=0; i<meshData.size(); i++) {
std::shared_ptr<IO::Mesh> mesh = meshData[i].mesh;
meshes_written.push_back( write_domain_silo(fid,filename,meshData[i],format) );
}
silo::close( fid );
return meshes_written;
#else
ERROR("Application built without silo support");
return std::vector<IO::MeshDatabase>();
#endif
}
/****************************************************
* Write the mesh data *
****************************************************/
void IO::writeData( const std::string& subdir, const std::vector<IO::MeshDataStruct>& meshData, MPI_Comm comm )
{
if ( global_IO_path.empty() )
IO::initialize( );
PROFILE_START("writeData");
int rank = comm_rank(comm);
// Create the output directory
char path[100];
sprintf(path,"vis%03i",timestep);
if ( rank == 0 )
mkdir(path,S_IRWXU|S_IRGRP);
std::string path = global_IO_path + "/" + subdir;
if ( rank == 0 ) {
mkdir(path.c_str(),S_IRWXU|S_IRGRP);
}
MPI_Barrier(comm);
// Write the mesh files
std::vector<IO::MeshDatabase> meshes_written;
if ( format == 1 ) {
if ( global_IO_format == Format::OLD ) {
// Write the original triangle format
meshes_written = writeMeshesOrigFormat( meshData, path );
} else if ( format == 2 ) {
// Write the new format
meshes_written = writeMeshesNewFormat( meshData, path, format );
} else if ( global_IO_format == Format::NEW ) {
// Write the new format (double precision)
meshes_written = writeMeshesNewFormat( meshData, path, 2 );
} else if ( global_IO_format == Format::SILO ) {
// Write silo
meshes_written = writeMeshesSilo( meshData, path, 4 );
} else {
ERROR("Unknown format");
}
@ -184,18 +351,29 @@ void IO::writeData( int timestep, const std::vector<IO::MeshDataStruct>& meshDat
if ( rank == 0 ) {
// Write the summary file for the current timestep
char filename[200];
sprintf(filename,"%s/LBM.summary",path);
sprintf(filename,"%s/LBM.summary",path.c_str());
write(meshes_written,filename);
// Add the timestep to the global summary file
FILE *fid = NULL;
if ( !global_summary_created ) {
fid = fopen("summary.LBM","wb");
global_summary_created = true;
} else {
fid = fopen("summary.LBM","ab");
// Write summary silo file if needed
#ifdef USE_SILO
if ( global_IO_format == Format::SILO ) {
sprintf(filename,"%s/summary.silo",path.c_str());
writeSiloSummary(meshes_written,filename);
}
#endif
// Add the timestep to the global summary file
if ( global_IO_format == Format::OLD || global_IO_format == Format::NEW ) {
auto filename = global_IO_path+"/summary.LBM";
FILE *fid = fopen(filename.c_str(),"ab");
fprintf(fid,"%s/\n",path.c_str());
fclose(fid);
} else if ( global_IO_format == Format::SILO ) {
auto filename = global_IO_path+"/LBM.visit";
FILE *fid = fopen(filename.c_str(),"ab");
fprintf(fid,"%s/summary.silo\n",path.c_str());
fclose(fid);
} else {
ERROR("Unknown format");
}
fprintf(fid,"%s\n",path);
fclose(fid);
}
PROFILE_STOP("writeData");
}

36
IO/Writer.h
View File

@ -12,18 +12,44 @@
namespace IO {
/*!
* @brief Initialize the writer
* @details This function initializes the writer to the given path. All subsequent
* writes will occur in this directory. If this is not called, then it will default
* to the current path.
* @param[in] path The path to use for writes
* @param[in] format The data format to use:
* old - Old mesh format (provided for backward compatibility, cannot write variables)
* new - New format, 1 file/process
* silo - Silo
* @param[in] append Append any existing data (default is false)
*/
void initialize( const std::string& path="", const std::string& format="new", bool append=false );
/*!
* @brief Write the data for the timestep
* @details This function writes the mesh and variable data provided for the current timestep
* @param[in] subdir The subdirectory to use for the timestep
* @param[in] meshData The data to write
* @param[in] comm The comm to use for writing (usually MPI_COMM_WORLD or a dup thereof)
*/
void writeData( const std::string& subdir, const std::vector<IO::MeshDataStruct>& meshData, MPI_Comm comm );
/*!
* @brief Write the data for the timestep
* @details This function writes the mesh and variable data provided for the current timestep
* @param[in] timestep The timestep iteration
* @param[in] meshData The data to write
* @param[in] format The data format to use:
* 1 - Old mesh format (provided for backward compatibility, cannot write variables)
* 2 - New format, 1 file/process, double precision
* 3 - New format, 1 file/process, single precision (not finished)
* @param[in] comm The comm to use for writing (usually MPI_COMM_WORLD or a dup thereof)
*/
void writeData( int timestep, const std::vector<IO::MeshDataStruct>& meshData, int format, MPI_Comm comm );
inline void writeData( int timestep, const std::vector<IO::MeshDataStruct>& meshData, MPI_Comm comm )
{
char subdir[100];
sprintf(subdir,"vis%03i",timestep);
writeData( subdir, meshData, comm );
}
} // IO namespace

214
IO/silo.cpp Normal file
View File

@ -0,0 +1,214 @@
#include "IO/silo.h"
#include "common/Utilities.h"
#include "common/MPI_Helpers.h"
#include "ProfilerApp.h"
#ifdef USE_SILO
#include <silo.h>
namespace silo {
/****************************************************
* Open/close a file *
****************************************************/
DBfile* open( const std::string& filename, FileMode mode )
{
DBfile *fid = nullptr;
if ( mode == CREATE ) {
fid = DBCreate( filename.c_str(), DB_CLOBBER, DB_LOCAL, nullptr, DB_HDF5 );
} else if ( mode == WRITE ) {
fid = DBOpen( filename.c_str(), DB_HDF5, DB_APPEND );
} else if ( mode == READ ) {
fid = DBOpen( filename.c_str(), DB_HDF5, DB_READ );
}
return fid;
}
void close( DBfile* fid )
{
DBClose( fid );
}
/****************************************************
* Write a uniform mesh to silo *
****************************************************/
template<int NDIM>
void writeUniformMesh( DBfile* fid, const std::string& meshname,
const std::array<double,2*NDIM>& range, const std::array<int,NDIM>& N )
{
PROFILE_START("writeUniformMesh",2);
int dims[NDIM];
for (size_t d=0; d<N.size(); d++)
dims[d] = N[d]+1;
float *x = nullptr;
if ( NDIM >= 1 ) {
x = new float[dims[0]];
for (int i=0; i<N[0]; i++)
x[i] = range[0] + i*(range[1]-range[0])/N[0];
x[N[0]] = range[1];
}
float *y = nullptr;
if ( NDIM >= 2 ) {
y = new float[dims[1]];
for (int i=0; i<N[1]; i++)
y[i] = range[2] + i*(range[3]-range[2])/N[1];
y[N[1]] = range[3];
}
float *z = nullptr;
if ( NDIM >= 3 ) {
z = new float[dims[2]];
for (int i=0; i<N[2]; i++)
z[i] = range[4] + i*(range[5]-range[4])/N[2];
z[N[2]] = range[5];
}
float *coords[] = { x, y, z };
int err = DBPutQuadmesh( fid, meshname.c_str(), nullptr, coords, dims, NDIM, DB_FLOAT, DB_COLLINEAR, nullptr );
ASSERT( err == 0 );
PROFILE_STOP("writeUniformMesh",2);
}
/****************************************************
* Write a vector/tensor quad variable *
****************************************************/
std::vector<std::string> getVarSuffix( int ndim, int nvars )
{
std::vector<std::string> suffix(nvars);
if ( nvars == 1 ) {
suffix[0] = "";
} else if ( nvars == ndim ) {
if ( ndim==2 ) {
suffix[0] = "_x";
suffix[1] = "_y";
} else if ( ndim==3 ) {
suffix[0] = "_x";
suffix[1] = "_y";
suffix[2] = "_z";
} else {
ERROR("Not finished");
}
} else if ( nvars == ndim*ndim ) {
if ( ndim==2 ) {
suffix[0] = "_xx";
suffix[1] = "_xy";
suffix[2] = "_yx";
suffix[3] = "_yy";
} else if ( ndim==3 ) {
suffix[0] = "_xx";
suffix[1] = "_xy";
suffix[2] = "_xz";
suffix[3] = "_yx";
suffix[4] = "_yy";
suffix[5] = "_yz";
suffix[6] = "_zx";
suffix[7] = "_zy";
suffix[8] = "_zz";
} else {
ERROR("Not finished");
}
} else {
for (int i=0; i<nvars; i++)
suffix[i] = "_" + std::to_string(i+1);
}
return suffix;
}
template<int NDIM>
void writeUniformMeshVariable( DBfile* fid, const std::string& meshname, const std::array<int,NDIM>& N,
const std::string& varname, const Array<double>& data, VariableType type )
{
for (int d=0; d<NDIM; d++)
ASSERT(N[d]==(int)data.size(d));
PROFILE_START("writeUniformMeshVariable",2);
int nvars=1, dims[NDIM]={1};
const double *vars[NDIM] = { nullptr };
if ( type == NodeVariable ) {
ERROR("Not finished");
} else if ( type == EdgeVariable ) {
ERROR("Not finished");
} else if ( type == SurfaceVariable ) {
ERROR("Not finished");
} else if ( type == VolumeVariable ) {
ASSERT( data.ndim()==NDIM || data.ndim()==NDIM+1 );
nvars = data.size(NDIM+1);
size_t N = data.length()/nvars;
for (int d=0; d<NDIM; d++)
dims[d] = data.size(d);
for (int i=0; i<nvars; i++)
vars[i] = &data(i*N);
} else {
ERROR("Invalid variable type");
}
auto suffix = getVarSuffix( NDIM, nvars );
std::vector<std::string> var_names(nvars);
for (int i=0; i<nvars; i++)
var_names[i] = varname + suffix[i];
std::vector<char*> varnames(nvars,nullptr);
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, DB_ZONECENT, nullptr );
ASSERT( err == 0 );
PROFILE_STOP("writeUniformMeshVariable",2);
}
/****************************************************
* Write a multimesh *
****************************************************/
void writeMultiMesh( DBfile* fid, const std::string& meshname,
const std::vector<std::string>& meshNames,
const std::vector<int>& meshTypes )
{
std::vector<char*> meshnames(meshNames.size());
for ( size_t i = 0; i < meshNames.size(); ++i )
meshnames[i] = (char *) meshNames[i].c_str();
std::string tree_name = meshname + "_tree";
DBoptlist *optList = DBMakeOptlist( 1 );
DBAddOption( optList, DBOPT_MRGTREE_NAME, (char *) tree_name.c_str() );
DBPutMultimesh( fid, meshname.c_str(), meshNames.size(), meshnames.data(), meshTypes.data(), nullptr );
DBFreeOptlist( optList );
}
/****************************************************
* Write a multivariable *
****************************************************/
void writeMultiVar( DBfile* fid, const std::string& varname,
const std::vector<std::string>& varNames,
const std::vector<int>& varTypes, int ndim, int nvar )
{
ASSERT(varNames.size()==varTypes.size());
auto suffix = getVarSuffix( ndim, nvar );
for (int i=0; i<nvar; i++) {
auto varname2 = varname + suffix[i];
auto varNames2 = varNames;
for ( auto& tmp : varNames2 )
tmp += suffix[i];
std::vector<char*> varnames(varNames.size(),nullptr);
for (size_t j=0; j<varNames.size(); j++)
varnames[j] = const_cast<char*>(varNames2[j].c_str());
DBPutMultivar( fid, varname2.c_str(), varNames.size(), varnames.data(), (int*) varTypes.data(), nullptr );
}
}
}; // silo namespace
// Explicit instantiations
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 );
#else
#endif

104
IO/silo.h Normal file
View File

@ -0,0 +1,104 @@
#ifndef NETCDF_READER
#define NETCDF_READER
#include <string>
#include <vector>
#include <array>
#include "common/Array.h"
#include "common/MPI_Helpers.h"
#include "common/Communication.h"
#ifdef USE_SILO
#include <silo.h>
#else
typedef int DBfile;
#endif
namespace silo {
enum FileMode { READ, WRITE, CREATE };
enum VariableType { NodeVariable=1, EdgeVariable=2, SurfaceVariable=2, VolumeVariable=3, NullVariable=0 };
/*!
* @brief Open silo file
* @detailed This function opens a silo file
* @param[in] filename File to open
* @param[in] mode Open the file for reading or writing
* @return This function returns a handle to the file
*/
DBfile* open( const std::string& filename, FileMode mode );
/*!
* @brief Close silo file
* @detailed This function closes a silo file
* @param[in] fid Handle to the open file
*/
void close( DBfile* fid );
/*!
* @brief Write a uniform grid
* @detailed This function writes a uniform grid to silo as a Quadmesh
* @return This function returns a handle to the file
* @param[in] meshname Mesh name
* @param[in] range Range of mesh { xmin, xmax, ymin, ymax, zmin, zmax }
* @param[in] N Number of cells in each direction
*/
template<int NDIM>
void writeUniformMesh( DBfile* fid, const std::string& meshname,
const std::array<double,2*NDIM>& range, const std::array<int,NDIM>& N );
/*!
* @brief Write a multimesh
* @detailed This function writes a multimesh to silo
* @return This function returns a handle to the file
* @param[in] meshname Mesh name
* @param[in] subMeshNames Names of the sub meshes in the form "filename:meshname"
* @param[in] subMeshTypes Type of each submesh
*/
void writeMultiMesh( DBfile* fid, const std::string& meshname,
const std::vector<std::string>& subMeshNames,
const std::vector<int>& subMeshTypes );
/*!
* @brief Write a uniform grid
* @detailed This function writes a uniform grid to silo as a Quadmesh
* @param[in] fid Handle to the open file
* @param[in] meshname Mesh name
* @param[in] N Number of cells in each direction
* @param[in] varname Variable name
* @param[in] data Variable data
* @param[in] type Variable type
*/
template<int NDIM>
void writeUniformMeshVariable( DBfile* fid, const std::string& meshname, const std::array<int,NDIM>& N,
const std::string& varname, const Array<double>& data, VariableType type );
/*!
* @brief Write a multivariable
* @detailed This function writes a multivariable to silo
* @return This function returns a handle to the file
* @param[in] fid Handle to the open file
* @param[in] varname Mesh name
* @param[in] subVarNames Names of the sub meshes in the form "filename:meshname"
* @param[in] subVarTypes Type of each submesh
* @param[in] ndim Dimension of variable (used to determine suffix)
* @param[in] nvar Number of subvariables (used to determine suffix)
*/
void writeMultiVar( DBfile* fid, const std::string& varname,
const std::vector<std::string>& subVarNames,
const std::vector<int>& subVarTypes, int ndim, int nvar );
}; // silo namespace
#endif

4
cmake/FindVisIt.cmake
View File

@ -48,7 +48,9 @@ MACRO( VISIT_PLUGIN SRC_DIR TARGET )
COMMAND ${VISIT_XML_CMAKE} -clobber ${TARGET}.xml
COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE}
-DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} -DCMAKE_C_FLAGS=${CMAKE_C_FLAGS}
-DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER} -DCMAKE_CXX_FLAGS=${CMAKE_CXX_FLAGS} .
-DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER} -DCMAKE_CXX_FLAGS=${CMAKE_CXX_FLAGS}
-DVISIT_NOLINK_MPI_WITH_LIBRARIES=TRUE
.
COMMAND make
WORKING_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}/${SRC_DIR}"
SOURCES ${SRC_DIR}

32
cmake/libraries.cmake
View File

@ -193,7 +193,7 @@ MACRO ( CONFIGURE_HDF5 )
${HDF5_HL_LIB}
${HDF5_LIB}
)
ADD_DEFINITIONS ( "-D USE_HDF5" )
ADD_DEFINITIONS ( -DUSE_HDF5 )
MESSAGE( "Using hdf5" )
MESSAGE( " ${HDF5_LIB}" )
ENDIF()
@ -227,13 +227,41 @@ MACRO( CONFIGURE_NETCDF )
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" )
ADD_DEFINITIONS ( -DUSE_NETCDF )
MESSAGE( "Using netcdf" )
MESSAGE( " ${NETCDF_LIBS}" )
ENDIF()
ENDMACRO()
# Macro to find and configure the silo libraries
MACRO ( CONFIGURE_SILO )
# Determine if we want to use silo
CHECK_ENABLE_FLAG( USE_EXT_SILO 0 )
IF ( USE_SILO )
SET( USE_HDF5 1 )
CONFIGURE_HDF5()
# Check if we specified the silo directory
IF ( SILO_DIRECTORY )
VERIFY_PATH ( ${SILO_DIRECTORY} )
INCLUDE_DIRECTORIES ( ${SILO_DIRECTORY}/include )
SET ( SILO_INCLUDE ${SILO_DIRECTORY}/include )
FIND_LIBRARY ( SILO_LIB NAMES siloh5 PATHS ${SILO_DIRECTORY}/lib NO_DEFAULT_PATH )
ELSE()
MESSAGE( "Default search for silo is not yet supported")
MESSAGE( "Use -D SILO_DIRECTORY=" FATAL_ERROR)
ENDIF()
SET ( SILO_LIBS
${SILO_LIB}
)
SET( EXTERNAL_LIBS ${EXTERNAL_LIBS} ${SILO_LIBS} ${HDF5_LIBS} )
ADD_DEFINITIONS ( -DUSE_SILO )
MESSAGE( "Using silo" )
MESSAGE( " ${SILO_LIB}" )
ENDIF ()
ENDMACRO()
# Macro to configure system-specific libraries and flags
MACRO( CONFIGURE_SYSTEM )
# First check/set the compile mode

2
common/TwoPhase.cpp
View File

@ -1086,7 +1086,7 @@ void TwoPhase::WriteSurfaces(int logcount)
meshData[1].mesh = ws_mesh;
meshData[2].meshName = "ns-tris";
meshData[2].mesh = ns_mesh;
IO::writeData( logcount, meshData, 2, Dm.Comm );
IO::writeData( logcount, meshData, Dm.Comm );
}

1
tests/CMakeLists.txt
View File

@ -19,6 +19,7 @@ ADD_LBPM_EXECUTABLE( ColorToBinary )
ADD_LBPM_EXECUTABLE( BlobAnalysis )
ADD_LBPM_EXECUTABLE( BlobIdentify )
ADD_LBPM_EXECUTABLE( BlobIdentifyParallel )
ADD_LBPM_EXECUTABLE( convertIO )
#ADD_LBPM_EXECUTABLE( BlobAnalyzeParallel )
CONFIGURE_FILE( ${CMAKE_CURRENT_SOURCE_DIR}/cylindertest ${CMAKE_CURRENT_BINARY_DIR}/cylindertest COPYONLY )

2
tests/ComponentLabel.cpp
View File

@ -476,7 +476,7 @@ int main(int argc, char **argv)
fillData.copy(Averages.Label_WP,LabelWPVar->data);
fillData.copy(Averages.Label_NWP,LabelNWPVar->data);
fillData.copy(Averages.PhaseID,PhaseIDVar->data);
IO::writeData( 0, meshData, 2, comm );
IO::writeData( 0, meshData, comm );
/*
FILE *NWP_FILE;
NWP_FILE = fopen("NWP.dat","wb");

8
tests/TestBlobIdentify.cpp
View File

@ -229,7 +229,7 @@ int main(int argc, char **argv)
fillData.copy(Phase,PhaseVar->data);
fillData.copy(SignDist,SignDistVar->data);
fillData.copy(GlobalBlobID,BlobIDVar->data);
IO::writeData( 0, meshData, 2, comm );
IO::writeData( 0, meshData, comm );
writeIDMap(ID_map_struct(nblobs),0,"lbpm_id_map.txt");
int save_it = 1;
@ -281,7 +281,7 @@ int main(int argc, char **argv)
fillData.copy(Phase,PhaseVar->data);
fillData.copy(SignDist,SignDistVar->data);
fillData.copy(GlobalBlobID,BlobIDVar->data);
IO::writeData( save_it, meshData, 2, comm );
IO::writeData( save_it, meshData, comm );
}
PROFILE_STOP("constant velocity test");
@ -305,7 +305,7 @@ int main(int argc, char **argv)
fillData.copy(Phase,PhaseVar->data);
fillData.copy(SignDist,SignDistVar->data);
fillData.copy(GlobalBlobID,BlobIDVar->data);
IO::writeData( save_it, meshData, 2, comm );
IO::writeData( save_it, meshData, comm );
save_it++;
id_max = nblobs-1;
for (int i=0; i<25; i++, save_it++) {
@ -404,7 +404,7 @@ int main(int argc, char **argv)
fillData.copy(Phase,PhaseVar->data);
fillData.copy(SignDist,SignDistVar->data);
fillData.copy(GlobalBlobID,BlobIDVar->data);
IO::writeData( save_it, meshData, 2, comm );
IO::writeData( save_it, meshData, comm );
}
PROFILE_STOP("moving bubble test");

2
tests/TestBlobIdentifyCorners.cpp
View File

@ -115,7 +115,7 @@ int main(int argc, char **argv)
fillData.copy(SignDist,SignDistVar->data);
fillData.copy(LocalBlobID,LocalBlobIDVar->data);
fillData.copy(GlobalBlobID,GlobalBlobIDVar->data);
IO::writeData( 0, meshData, 2 );
IO::writeData( 0, meshData );
int N_errors = 0;

2
tests/TestBubble.cpp
View File

@ -2187,7 +2187,7 @@ int main(int argc, char **argv)
}
meshData[k].vars.push_back(dist);
}
IO::writeData( bubbleCount, meshData, 2, comm );
IO::writeData( bubbleCount, meshData, comm );
#else
fclose(WN_TRIS);
#endif

317
tests/TestWriter.cpp
View File

@ -31,7 +31,181 @@ 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 )
{
int rank, nprocs;
MPI_Comm comm = MPI_COMM_WORLD;
MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs);
MPI_Barrier(comm);
// Write the data
PROFILE_START(format+"-write");
IO::initialize( "test_"+format, format, false );
IO::writeData( 0, meshData, comm );
IO::writeData( 3, meshData, comm );
MPI_Barrier(comm);
PROFILE_STOP(format+"-write");
// Get the summary name for reading
std::string summary_name;
if ( format=="old" || format=="new" )
summary_name = "test_" + format + "/summary.LBM";
else if ( format == "silo" )
summary_name = "test_" + format + "/LBM.visit";
else
ERROR("Unknown format");
// Get direct access to the meshes used to test the reader
const auto pointmesh = dynamic_cast<IO::PointList*>( meshData[0].mesh.get() );
const auto trimesh = dynamic_cast<IO::TriMesh*>( meshData[1].mesh.get() );
const auto trilist = dynamic_cast<IO::TriList*>( meshData[2].mesh.get() );
const auto domain = dynamic_cast<IO::DomainMesh*>( meshData[3].mesh.get() );
const size_t N_tri = trimesh->A.size();
// Get a list of the timesteps
PROFILE_START(format+"-read-timesteps");
std::vector<std::string> timesteps = IO::readTimesteps(summary_name);
PROFILE_STOP(format+"-read-timesteps");
if ( timesteps.size()==2 )
ut.passes("Corrent number of timesteps");
else
ut.failure("Incorrent number of timesteps");
// Check the mesh lists
for (size_t i=0; i<timesteps.size(); i++) {
// Load the list of meshes and check its size
PROFILE_START(format+"-read-getMeshList");
std::vector<IO::MeshDatabase> list = IO::getMeshList(".",timesteps[i]);
PROFILE_STOP(format+"-read-getMeshList");
if ( list.size()==meshData.size() )
ut.passes("Corrent number of meshes found");
else
ut.failure("Incorrent number of meshes found");
// Check the number of domains for each mesh
bool pass = true;
for (size_t j=0; j<list.size(); j++)
pass = pass && (int)list[j].domains.size()==nprocs;
if ( pass ) {
ut.passes("Corrent number of domains for mesh");
} else {
ut.failure("Incorrent number of domains for mesh");
continue;
}
// For each domain, load the mesh and check its data
for (size_t j=0; j<list.size(); j++) {
for (size_t k=0; k<list[j].domains.size(); k++) {
PROFILE_START(format+"-read-getMesh");
std::shared_ptr<IO::Mesh> mesh = IO::getMesh(".",timesteps[i],list[j],k);
PROFILE_STOP(format+"-read-getMesh");
if ( mesh.get()==NULL ) {
printf("Failed to load %s\n",list[j].name.c_str());
pass = false;
break;
}
if ( list[j].name=="pointmesh" ) {
// Check the pointmesh
std::shared_ptr<IO::PointList> pmesh = IO::getPointList(mesh);
if ( pmesh.get()==NULL ) {
pass = false;
break;
}
if ( pmesh->points.size() != pointmesh->points.size() ) {
pass = false;
break;
}
}
if ( list[j].name=="trimesh" || list[j].name=="trilist" ) {
// Check the trimesh/trilist
std::shared_ptr<IO::TriMesh> mesh1 = IO::getTriMesh(mesh);
std::shared_ptr<IO::TriList> mesh2 = IO::getTriList(mesh);
if ( mesh1.get()==NULL || mesh2.get()==NULL ) {
pass = false;
break;
}
if ( mesh1->A.size()!=N_tri || mesh1->B.size()!=N_tri || mesh1->C.size()!=N_tri ||
mesh2->A.size()!=N_tri || mesh2->B.size()!=N_tri || mesh2->C.size()!=N_tri )
{
pass = false;
break;
}
const std::vector<Point>& P1 = mesh1->vertices->points;
const std::vector<int>& A1 = mesh1->A;
const std::vector<int>& B1 = mesh1->B;
const std::vector<int>& C1 = mesh1->C;
const std::vector<Point>& A2 = mesh2->A;
const std::vector<Point>& B2 = mesh2->B;
const std::vector<Point>& C2 = mesh2->C;
const std::vector<Point>& A = trilist->A;
const std::vector<Point>& B = trilist->B;
const std::vector<Point>& C = trilist->C;
for (size_t i=0; i<N_tri; i++) {
if ( !approx_equal(P1[A1[i]],A[i]) || !approx_equal(P1[B1[i]],B[i]) || !approx_equal(P1[C1[i]],C[i]) )
pass = false;
if ( !approx_equal(A2[i],A[i]) || !approx_equal(B2[i],B[i]) || !approx_equal(C2[i],C[i]) )
pass = false;
}
}
if ( list[j].name=="domain" && format!="old" ) {
// Check the domain mesh
const IO::DomainMesh& mesh1 = *std::dynamic_pointer_cast<IO::DomainMesh>(mesh);
if ( mesh1.nprocx!=domain->nprocx || mesh1.nprocy!=domain->nprocy || mesh1.nprocz!=domain->nprocz )
pass = false;
if ( mesh1.nx!=domain->nx || mesh1.ny!=domain->ny || mesh1.nz!=domain->nz )
pass = false;
if ( mesh1.Lx!=domain->Lx || mesh1.Ly!=domain->Ly || mesh1.Lz!=domain->Lz )
pass = false;
}
}
}
if ( pass ) {
ut.passes("Meshes loaded correctly");
} else {
ut.failure("Meshes did not load correctly");
continue;
}
// For each domain, load the variables and check their data
if ( format=="old" )
continue; // Old format does not support variables
for (size_t j=0; j<list.size(); j++) {
const IO::MeshDataStruct* mesh0 = NULL;
for (size_t k=0; k<meshData.size(); k++) {
if ( meshData[k].meshName == list[j].name ) {
mesh0 = &meshData[k];
break;
}
}
for (size_t v=0; v<list[j].variables.size(); v++) {
for (size_t k=0; k<list[j].domains.size(); k++) {
PROFILE_START(format+"-read-getVariable");
std::shared_ptr<const IO::Variable> variable =
IO::getVariable(".",timesteps[i],list[j],k,list[j].variables[v].name);
PROFILE_STOP(format+"-read-getVariable");
const IO::Variable& var1 = *mesh0->vars[v];
const IO::Variable& var2 = *variable;
pass = var1.name == var2.name;
pass = pass && var1.dim == var2.dim;
pass = pass && var1.type == var2.type;
pass = pass && var1.data.length() == var2.data.length();
if ( pass ) {
for (size_t m=0; m<var1.data.length(); m++)
pass = pass && approx_equal(var1.data(m),var2.data(m));
}
if ( pass ) {
ut.passes("Variables match");
} else {
ut.failure("Variables did not match");
break;
}
}
}
}
}
}
// Main
int main(int argc, char **argv)
{
int rank,nprocs;
@ -128,147 +302,14 @@ int main(int argc, char **argv)
meshData[3].mesh = domain;
meshData[3].vars.push_back(domain_var);
// Write the data
std::vector<int> format(2,0);
format[0] = 2;
format[1] = 1;
IO::writeData( 0, meshData, format[0], comm );
IO::writeData( 3, meshData, format[1], comm );
MPI_Barrier(comm);
// Get a list of the timesteps
std::vector<std::string> timesteps = IO::readTimesteps("summary.LBM");
if ( timesteps.size()==2 )
ut.passes("Corrent number of timesteps");
else
ut.failure("Incorrent number of timesteps");
// Check the mesh lists
for (size_t i=0; i<timesteps.size(); i++) {
// Load the list of meshes and check its size
std::vector<IO::MeshDatabase> list = IO::getMeshList(".",timesteps[i]);
if ( list.size()==meshData.size() )
ut.passes("Corrent number of meshes found");
else
ut.failure("Incorrent number of meshes found");
// Check the number of domains for each mesh
bool pass = true;
for (size_t j=0; j<list.size(); j++)
pass = pass && (int)list[j].domains.size()==nprocs;
if ( pass ) {
ut.passes("Corrent number of domains for mesh");
} else {
ut.failure("Incorrent number of domains for mesh");
continue;
}
// For each domain, load the mesh and check its data
for (size_t j=0; j<list.size(); j++) {
for (size_t k=0; k<list[j].domains.size(); k++) {
std::shared_ptr<IO::Mesh> mesh = IO::getMesh(".",timesteps[i],list[j],k);
if ( mesh.get()==NULL ) {
printf("Failed to load %s\n",list[j].name.c_str());
pass = false;
break;
}
if ( list[j].name=="pointmesh" ) {
// Check the pointmesh
std::shared_ptr<IO::PointList> pmesh = IO::getPointList(mesh);
if ( pmesh.get()==NULL ) {
pass = false;
break;
}
if ( pmesh->points.size()!=N_points ) {
pass = false;
break;
}
}
if ( list[j].name=="trimesh" || list[j].name=="trilist" ) {
// Check the trimesh/trilist
std::shared_ptr<IO::TriMesh> mesh1 = IO::getTriMesh(mesh);
std::shared_ptr<IO::TriList> mesh2 = IO::getTriList(mesh);
if ( mesh1.get()==NULL || mesh2.get()==NULL ) {
pass = false;
break;
}
if ( mesh1->A.size()!=N_tri || mesh1->B.size()!=N_tri || mesh1->C.size()!=N_tri ||
mesh2->A.size()!=N_tri || mesh2->B.size()!=N_tri || mesh2->C.size()!=N_tri )
{
pass = false;
break;
}
const std::vector<Point>& P1 = mesh1->vertices->points;
const std::vector<int>& A1 = mesh1->A;
const std::vector<int>& B1 = mesh1->B;
const std::vector<int>& C1 = mesh1->C;
const std::vector<Point>& A2 = mesh2->A;
const std::vector<Point>& B2 = mesh2->B;
const std::vector<Point>& C2 = mesh2->C;
const std::vector<Point>& A = trilist->A;
const std::vector<Point>& B = trilist->B;
const std::vector<Point>& C = trilist->C;
for (size_t i=0; i<N_tri; i++) {
if ( !approx_equal(P1[A1[i]],A[i]) || !approx_equal(P1[B1[i]],B[i]) || !approx_equal(P1[C1[i]],C[i]) )
pass = false;
if ( !approx_equal(A2[i],A[i]) || !approx_equal(B2[i],B[i]) || !approx_equal(C2[i],C[i]) )
pass = false;
}
}
if ( list[j].name=="domain" && format[i]>2 ) {
// Check the domain mesh
const IO::DomainMesh& mesh1 = *std::dynamic_pointer_cast<IO::DomainMesh>(mesh);
if ( mesh1.nprocx!=domain->nprocx || mesh1.nprocy!=domain->nprocy || mesh1.nprocz!=domain->nprocz )
pass = false;
if ( mesh1.nx!=domain->nx || mesh1.ny!=domain->ny || mesh1.nz!=domain->nz )
pass = false;
if ( mesh1.Lx!=domain->Lx || mesh1.Ly!=domain->Ly || mesh1.Lz!=domain->Lz )
pass = false;
}
}
}
if ( pass ) {
ut.passes("Meshes loaded correctly");
} else {
ut.failure("Meshes did not load correctly");
continue;
}
// For each domain, load the variables and check their data
if ( format[i]==1 )
continue; // Format 1 does not support variables
for (size_t j=0; j<list.size(); j++) {
const IO::MeshDataStruct* mesh0 = NULL;
for (size_t k=0; k<meshData.size(); k++) {
if ( meshData[k].meshName == list[j].name ) {
mesh0 = &meshData[k];
break;
}
}
for (size_t v=0; v<list[j].variables.size(); v++) {
for (size_t k=0; k<list[j].domains.size(); k++) {
std::shared_ptr<const IO::Variable> variable =
IO::getVariable(".",timesteps[i],list[j],k,list[j].variables[v].name);
const IO::Variable& var1 = *mesh0->vars[v];
const IO::Variable& var2 = *variable;
pass = var1.name == var2.name;
pass = pass && var1.dim == var2.dim;
pass = pass && var1.type == var2.type;
pass = pass && var1.data.length() == var2.data.length();
if ( pass ) {
for (size_t m=0; m<var1.data.length(); m++)
pass = pass && approx_equal(var1.data(m),var2.data(m));
}
if ( pass ) {
ut.passes("Variables match");
} else {
ut.failure("Variables did not match");
break;
}
}
}
}
}
// Run the tests
testWriter( "old", meshData, ut );
testWriter( "new", meshData, ut );
//testWriter( "silo", meshData, ut );
// Finished
ut.report();
PROFILE_SAVE("TestWriter",true);
int N_errors = ut.NumFailGlobal();
MPI_Barrier(comm);
MPI_Finalize();

90
tests/convertIO.cpp Normal file
View File

@ -0,0 +1,90 @@
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <exception>
#include <stdexcept>
#include <fstream>
#include "common/MPI_Helpers.h"
#include "common/Communication.h"
#include "common/Utilities.h"
#include "IO/Mesh.h"
#include "IO/Reader.h"
#include "IO/Writer.h"
#include "ProfilerApp.h"
int main(int argc, char **argv)
{
// Initialize MPI
int rank,nprocs;
MPI_Init(&argc,&argv);
MPI_Comm comm = MPI_COMM_WORLD;
MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs);
Utilities::setErrorHandlers();
PROFILE_ENABLE(2);
PROFILE_ENABLE_TRACE();
PROFILE_START("Main");
{ // Limit scope
// Get inputs
if ( argc != 3 ) {
std::cerr << "Error calling convertIO:\n";
std::cerr << " convertIO input_file format\n";
return -1;
}
std::string filename = argv[1];
std::string format = argv[2];
std::string path = IO::getPath( filename );
// Read the timesteps
auto timesteps = IO::readTimesteps( filename );
// Loop through the timesteps, reading/writing the data
IO::initialize( "", format, false );
for ( auto timestep : timesteps ) {
// Read the list of MeshDatabase
auto databases = IO::getMeshList( path, timestep );
// Build the MeshDataStruct
std::vector<IO::MeshDataStruct> meshData(databases.size());
// Loop through the database
int i = 0;
PROFILE_START("Read");
for ( const auto& database : databases ) {
// Read the appropriate mesh domain
ASSERT( (int) database.domains.size() == nprocs );
meshData[i].meshName = database.name;
meshData[i].mesh = IO::getMesh( path, timestep, database, rank );
// Read the variables
for ( auto var : database.variables ) {
auto varData = IO::getVariable( path, timestep, database, rank, var.name );
IO::reformatVariable( *meshData[i].mesh, *varData );
meshData[i].vars.push_back( varData );
}
i++;
}
MPI_Barrier(comm);
PROFILE_STOP("Read");
// Save the mesh data to a new file
PROFILE_START("Write");
IO::writeData( timestep, meshData, MPI_COMM_WORLD );
MPI_Barrier(comm);
PROFILE_STOP("Write");
}
} // Limit scope
PROFILE_STOP("Main");
PROFILE_SAVE("convertData",true);
MPI_Barrier(comm);
MPI_Finalize();
return 0;
}

2
tests/lbpm_color_simulator.h
View File

@ -164,7 +164,7 @@ public:
fillData.copy(Averages.Press,PressData);
fillData.copy(Averages.SDs,SignData);
fillData.copy(Averages.Label_NWP,BlobData);
IO::writeData( timestep, visData, 2, newcomm );
IO::writeData( timestep, visData, newcomm );
PROFILE_STOP("Save Vis",1);
ThreadPool::WorkItem::d_state = 2; // Change state to finished
};

4
tests/lbpm_morphopen_pp.cpp
View File

@ -206,7 +206,7 @@ int main(int argc, char **argv)
// Compute min and max poresize
if (rank==0) printf(" computing local minimum and maximum... \n");
MinPoreSize=MaxPoreSize=PoreSize[0];
for (int idx=0; idx<PoreSize.size(); idx++){
for (size_t idx=0; idx<PoreSize.size(); idx++){
if (PoreSize[idx] < MinPoreSize) MinPoreSize=PoreSize[idx];
if (PoreSize[idx] > MaxPoreSize) MaxPoreSize=PoreSize[idx];
}
@ -220,7 +220,7 @@ int main(int argc, char **argv)
// Generate histogram counts
if (rank==0) printf(" generating local bin counts... \n");
BinWidth=(MaxPoreSize-MinPoreSize)/double(NumBins);
for (int idx=0; idx<PoreSize.size(); idx++){
for (size_t idx=0; idx<PoreSize.size(); idx++){
double value = PoreSize[idx];
int myBin = 0;
while (MinPoreSize+myBin*BinWidth < value) myBin++;

2
tests/lbpm_uCT_maskfilter.cpp
View File

@ -392,7 +392,7 @@ int main(int argc, char **argv)
#endif
// Write visulization data
IO::writeData( 0, meshData, 2, comm );
IO::writeData( 0, meshData, comm );
if (rank==0) printf("Finished. \n");

2
tests/lbpm_uCT_pp.cpp
View File

@ -354,7 +354,7 @@ int main(int argc, char **argv)
#endif
// Write visulization data
IO::writeData( 0, meshData, 2, comm );
IO::writeData( 0, meshData, comm );
if (rank==0) printf("Finished. \n");
/* for (k=0;k<nz;k++){