Commit eeccd0b6 authored by Burlen Loring's avatar Burlen Loring
Browse files

posthocIO analysis adapter

I added an adapter to perform I/O for conventional
post process analysis. It currently supports MPI-I/O
collective buffering.
parent 05bc16a3
<sensei>
<analysis type="PosthocIO" array="data" association="cell" enabled="0" period="2" />
<analysis type="histogram" array="data" association="cell" bins="10" />
<analysis type="autocorrelation" array="data" association="cell" window="10" k-max="3" />
<analysis type="catalyst" pipeline="slice" array="data" association="cell" />
......
......@@ -174,7 +174,7 @@ void analysis_final(size_t k_max, size_t nblocks)
std::cout << ' ' << result[i];
std::cout << std::endl;
}
master->foreach<Autocorrelation>([](Autocorrelation* b, const diy::Master::ProxyWithLink& cp, void*)
master->foreach<Autocorrelation>([](Autocorrelation*, const diy::Master::ProxyWithLink& cp, void*)
{
cp.collectives()->clear();
});
......@@ -183,10 +183,10 @@ void analysis_final(size_t k_max, size_t nblocks)
diy::ContiguousAssigner assigner(master->communicator().size(), nblocks); // NB: this is coupled to main(...) in oscillator.cpp
diy::RegularMergePartners partners(3, nblocks, 2, true);
diy::reduce(*master, assigner, partners,
[k_max](void* b_, const diy::ReduceProxy& rp, const diy::RegularMergePartners& partners)
[k_max](void* b_, const diy::ReduceProxy& rp, const diy::RegularMergePartners&)
{
Autocorrelation* b = static_cast<Autocorrelation*>(b_);
unsigned round = rp.round(); // current round number
//unsigned round = rp.round(); // current round number
using MaxHeapVector = std::vector<std::vector<std::tuple<float,Vertex>>>;
using Compare = std::greater<std::tuple<float,Vertex>>;
......@@ -211,7 +211,7 @@ void analysis_final(size_t k_max, size_t nblocks)
});
} else
{
for (size_t i = 0; i < rp.in_link().size(); ++i)
for (long i = 0; i < rp.in_link().size(); ++i)
{
MaxHeapVector in_heaps;
rp.dequeue(rp.in_link().target(i).gid, in_heaps);
......
......@@ -24,9 +24,11 @@ void initialize(MPI_Comm world,
int* to_x, int* to_y, int* to_z,
const std::string& config_file)
{
(void)window;
GlobalDataAdaptor = vtkSmartPointer<oscillators::DataAdaptor>::New();
GlobalDataAdaptor->Initialize(nblocks);
GlobalDataAdaptor->SetDataTimeStep(-1);
for (size_t cc=0; cc < n_local_blocks; ++cc)
{
GlobalDataAdaptor->SetBlockExtent(gid[cc],
......@@ -35,6 +37,9 @@ void initialize(MPI_Comm world,
from_z[cc], to_z[cc]);
}
int dext[6] = {0, domain_shape_x, 0, domain_shape_y, 0, domain_shape_z};
GlobalDataAdaptor->SetDataExtent(dext);
GlobalAnalysisAdaptor = vtkSmartPointer<sensei::ConfigurableAnalysis>::New();
GlobalAnalysisAdaptor->Initialize(world, config_file);
}
......@@ -57,6 +62,8 @@ void analyze(float time)
//-----------------------------------------------------------------------------
void finalize(size_t k_max, size_t nblocks)
{
(void)k_max;
(void)nblocks;
GlobalAnalysisAdaptor = NULL;
GlobalDataAdaptor = NULL;
}
......
#include "dataadaptor.h"
#include <vtkInformation.h>
#include <vtkFloatArray.h>
#include <vtkImageData.h>
#include <vtkMultiBlockDataSet.h>
......@@ -12,13 +13,13 @@
namespace oscillators
{
class DataAdaptor::DInternals
struct DataAdaptor::DInternals
{
public:
std::vector<diy::DiscreteBounds> CellExtents;
std::vector<float*> Data;
vtkSmartPointer<vtkMultiBlockDataSet> Mesh;
std::vector<vtkSmartPointer<vtkImageData> > BlockMesh;
std::vector<int> DataExtent;
};
inline bool areBoundsValid(const diy::DiscreteBounds& bds)
......@@ -74,6 +75,15 @@ void DataAdaptor::SetBlockExtent(int gid,
internals.CellExtents[gid].max[2] = zmax;
}
//-----------------------------------------------------------------------------
void DataAdaptor::SetDataExtent(int ext[6])
{
// TODO -- this key holds a int**, it should copy the data
this->Internals->DataExtent.assign(ext, ext+6);
this->GetInformation()->Set(vtkDataObject::DATA_EXTENT(),
this->Internals->DataExtent.data(), 6);
}
//-----------------------------------------------------------------------------
void DataAdaptor::SetBlockData(int gid, float* data)
{
......@@ -94,6 +104,8 @@ vtkDataObject* DataAdaptor::GetMesh(bool vtkNotUsed(structure_only))
internals.Mesh->SetBlock(static_cast<unsigned int>(cc), this->GetBlockMesh(cc));
}
}
this->AddArray(internals.Mesh,
vtkDataObject::FIELD_ASSOCIATION_CELLS, "data");
return internals.Mesh;
}
......@@ -117,13 +129,11 @@ vtkDataObject* DataAdaptor::GetBlockMesh(int gid)
//-----------------------------------------------------------------------------
bool DataAdaptor::AddArray(vtkDataObject* mesh, int association, const char* arrayname)
{
#ifndef NDEBUG
if (association != vtkDataObject::FIELD_ASSOCIATION_CELLS ||
arrayname == NULL ||
strcmp(arrayname, "data") != 0)
{
arrayname == NULL || strcmp(arrayname, "data") != 0)
return false;
}
#endif
bool retVal = false;
DInternals& internals = (*this->Internals);
vtkMultiBlockDataSet* md = vtkMultiBlockDataSet::SafeDownCast(mesh);
......
......@@ -21,13 +21,16 @@ public:
void SetBlockExtent(int gid,
int xmin, int xmax, int ymin, int ymax, int zmin, int zmax);
/// @brief Set the extent of the dataset.
void SetDataExtent(int ext[6]);
/// Set data for a specific block.
void SetBlockData(int gid, float* data);
virtual vtkDataObject* GetMesh(bool structure_only=false);
virtual bool AddArray(vtkDataObject* mesh, int association, const char* arrayname);
virtual unsigned int GetNumberOfArrays(int association) { return 1; }
virtual const char* GetArrayName(int association, unsigned int index)
virtual unsigned int GetNumberOfArrays(int) { return 1; }
virtual const char* GetArrayName(int, unsigned int index)
{ return index==0? "data" : NULL; }
virtual void ReleaseData();
......
......@@ -138,7 +138,7 @@ int main(int argc, char** argv)
// decompose the domain
diy::decompose(3, world.rank(), domain, assigner,
[&](int gid, const Bounds& core, const Bounds& bounds, const Bounds& domain, const Link& link)
[&](int gid, const Bounds&, const Bounds& bounds, const Bounds& domain, const Link& link)
{
auto b = new Block(gid, bounds, { domain.max[0] + 1, domain.max[1] + 1, domain.max[2] + 1 }, oscillators);
master.add(gid, b, new Link(link));
......
......@@ -28,6 +28,12 @@ void DataAdaptor::Initialize(
int tot_blocks_x, int tot_blocks_y, int tot_blocks_z,
int block_id_x, int block_id_y, int block_id_z)
{
(void)tot_blocks_x;
(void)tot_blocks_y;
(void)tot_blocks_z;
(void)block_id_x;
(void)block_id_y;
(void)block_id_z;
// we only really need to save the local extents for our current example. So
// we'll just save that.
this->CellExtent[0] = start_extents_x;
......
......@@ -296,7 +296,7 @@ void Autocorrelation::PrintResults(size_t k_max)
}
internals.Master->foreach<AutocorrelationImpl>(
[](AutocorrelationImpl* b, const diy::Master::ProxyWithLink& cp, void*)
[](AutocorrelationImpl*, const diy::Master::ProxyWithLink& cp, void*)
{
cp.collectives()->clear();
});
......@@ -305,10 +305,10 @@ void Autocorrelation::PrintResults(size_t k_max)
diy::ContiguousAssigner assigner(internals.Master->communicator().size(), nblocks); // NB: this is coupled to main(...) in oscillator.cpp
diy::RegularMergePartners partners(3, nblocks, 2, true);
diy::reduce(*internals.Master, assigner, partners,
[k_max](void* b_, const diy::ReduceProxy& rp, const diy::RegularMergePartners& partners)
[k_max](void* b_, const diy::ReduceProxy& rp, const diy::RegularMergePartners&)
{
AutocorrelationImpl* b = static_cast<AutocorrelationImpl*>(b_);
unsigned round = rp.round(); // current round number
AutocorrelationImpl* b = static_cast<AutocorrelationImpl*>(b_);
//unsigned round = rp.round(); // current round number
using MaxHeapVector = std::vector<std::vector<std::tuple<float,Vertex>>>;
using Compare = std::greater<std::tuple<float,Vertex>>;
......@@ -333,7 +333,7 @@ void Autocorrelation::PrintResults(size_t k_max)
});
} else
{
for (size_t i = 0; i < rp.in_link().size(); ++i)
for (long i = 0; i < rp.in_link().size(); ++i)
{
MaxHeapVector in_heaps;
rp.dequeue(rp.in_link().target(i).gid, in_heaps);
......
......@@ -8,7 +8,9 @@ set (sources
ConfigurableAnalysis.h
DataAdaptor.cxx
DataAdaptor.h
)
PosthocIO.cxx
PosthocIO.h
)
if(VTK_HAS_GENERIC_ARRAYS)
list(APPEND sources
......@@ -67,11 +69,12 @@ if(ENABLE_ADIOS)
endif()
target_link_libraries(sensei
PUBLIC ${VTK_LIBRARIES}
PUBLIC ${VTK_LIBRARIES}
mpi
PRIVATE diy
grid
pugixml
ArrayIO
)
#------------------------------------------------------------------------------
......
......@@ -6,6 +6,7 @@
#include <vtkDataObject.h>
#include "Autocorrelation.h"
#include "PosthocIO.h"
#ifdef ENABLE_HISTOGRAM
# include "Histogram.h"
#endif
......@@ -19,6 +20,11 @@
#include <vector>
#include <pugixml.hpp>
#define sensei_error(_arg) \
cerr << "ERROR: " << __FILE__ " : " << __LINE__ << std::endl \
<< "" _arg << std::endl;
namespace sensei
{
......@@ -35,7 +41,7 @@ class ConfigurableAnalysis::vtkInternals
{
return vtkDataObject::FIELD_ASSOCIATION_CELLS;
}
cout << "Invalid association type '" << association.c_str() << "'. Assuming 'point'" << endl;
sensei_error(<< "Invalid association type '" << association.c_str() << "'. Assuming 'point'");
return vtkDataObject::FIELD_ASSOCIATION_POINTS;
}
......@@ -57,7 +63,7 @@ public:
}
else
{
cerr << "'histogram' missing required attribute 'array'. Skipping." << endl;
sensei_error(<< "'histogram' missing required attribute 'array'. Skipping.");
}
}
#endif
......@@ -113,9 +119,65 @@ public:
node.attribute("k-max")? node.attribute("k-max").as_int() : 3);
this->Analyses.push_back(adaptor.GetPointer());
}
void AddPosthocIO(MPI_Comm comm, pugi::xml_node node)
{
if (!node.attribute("enabled") || !node.attribute("enabled").as_int())
{
cerr << "Skipping 'PosthocIO'." << endl;
return;
}
if (!node.attribute("array"))
{
sensei_error(<< "need at least one array");
return;
}
std::string arrayName = node.attribute("array").value();
std::vector<std::string> cellArrays;
std::vector<std::string> pointArrays;
pugi::xml_attribute assoc_att = node.attribute("association");
if (assoc_att && (std::string(assoc_att.value()) == "cell"))
cellArrays.push_back(arrayName);
else
pointArrays.push_back(arrayName);
std::string outputDir = "./";
if (node.attribute("output_dir"))
outputDir = node.attribute("output_dir").value();
std::string fileBase = "PosthocIO";
if (node.attribute("file_base"))
fileBase = node.attribute("file_base").value();
std::string blockExt = "block";
if (node.attribute("block_ext"))
blockExt = node.attribute("block_ext").value();
int mode = PosthocIO::MpiIO;
if (node.attribute("mode"))
mode = node.attribute("mode").as_int();
int period = 1;
if (node.attribute("period"))
period = node.attribute("period").as_int();
PosthocIO *adapter = PosthocIO::New();
adapter->Initialize(comm, outputDir, fileBase,
blockExt, cellArrays, pointArrays, mode,
period);
this->Analyses.push_back(adapter);
adapter->Delete();
}
};
//----------------------------------------------------------------------------
vtkStandardNewMacro(ConfigurableAnalysis);
//----------------------------------------------------------------------------
ConfigurableAnalysis::ConfigurableAnalysis()
: Internals(new ConfigurableAnalysis::vtkInternals())
......@@ -170,8 +232,16 @@ bool ConfigurableAnalysis::Initialize(MPI_Comm world, const std::string& filenam
this->Internals->AddAutoCorrelation(world, analysis);
continue;
}
cerr << "Skipping '" << type.c_str() << "'." << endl;
if (type == "PosthocIO")
{
this->Internals->AddPosthocIO(world, analysis);
continue;
}
std::cerr << "Skipping '" << type.c_str() << "'." << std::endl;
}
return true;
}
//----------------------------------------------------------------------------
......
#include "PosthocIO.h"
#include <vtkCompositeDataIterator.h>
#include <vtkCompositeDataSet.h>
#include <vtkDataArray.h>
#include <vtkDataObject.h>
#include <vtkDataSetAttributes.h>
#include <vtkObjectFactory.h>
#include <vtkSmartPointer.h>
#include <vtkImageData.h>
#include <vtkDataArrayTemplate.txx>
#include <vtkCellData.h>
#include <vtkPointData.h>
#include "DataAdaptor.h"
#include <algorithm>
#include <sstream>
#include <fstream>
#include <ArrayIO.h>
//#define PosthocIO_DEBUG
#define posthocIO_error(_arg) \
cerr << "ERROR: " << __FILE__ " : " << __LINE__ << std::endl \
<< "" _arg << std::endl;
#define posthocIO_status(_cond, _arg) \
if (_cond) \
{ \
cerr << "STATUS: " << __FILE__ " : " << __LINE__ << std::endl \
<< "" _arg << std::endl; \
}
namespace {
// **************************************************************************
void getWholePointExtents(vtkInformation *info, int *wholePointExtent)
{
memcpy(wholePointExtent, info->Get(vtkDataObject::DATA_EXTENT()),
6*sizeof(int));
}
// **************************************************************************
void getWholeCellExtents(vtkInformation *info, int *wholeCellExtent)
{
int *wholePointExtent = info->Get(vtkDataObject::DATA_EXTENT());
for (int i = 0; i < 6; ++i)
wholeCellExtent[i] = wholePointExtent[i] - (i%2 ? 1 : 0);
}
// **************************************************************************
void getLocalCellExtents(vtkImageData *id, int *localCellExtent)
{
id->GetExtent(localCellExtent);
localCellExtent[1] -= 1;
localCellExtent[3] -= 1;
localCellExtent[5] -= 1;
}
/*/ **************************************************************************
void getValidCellExtents(vtkImageData *id, int *validCellExtent)
{ getLocalCellExtents(id, validCellExtent); } */
// **************************************************************************
void getLocalPointExtents(vtkImageData *id, int *localPointExtent)
{
id->GetExtent(localPointExtent);
}
// **************************************************************************
void getValidPointExtents(vtkImageData *id,
int *wholePointExtent, int *validPointExtent)
{
// deal with coincident block faces
id->GetExtent(validPointExtent);
validPointExtent[1] = (validPointExtent[1] == wholePointExtent[1]
? validPointExtent[1] : validPointExtent[1] - 1);
validPointExtent[3] = (validPointExtent[3] == wholePointExtent[3]
? validPointExtent[3] : validPointExtent[3] - 1);
validPointExtent[5] = (validPointExtent[5] == wholePointExtent[5]
? validPointExtent[5] : validPointExtent[5] - 1);
}
// ****************************************************************************
int write(MPI_File file, MPI_Info hints,
int domain[6], int decomp[6], int valid[6], vtkDataArray *da,
bool useCollectives)
{
switch (da->GetDataType())
{
vtkTemplateMacro(
vtkDataArrayTemplate<VTK_TT> *ta =
static_cast<vtkDataArrayTemplate<VTK_TT>*>(da);
if ((useCollectives && arrayIO::write_all(file, hints,
domain, decomp, valid, ta->GetPointer(0))) ||
arrayIO::write(file, hints, domain, decomp,
valid, ta->GetPointer(0)))
{
posthocIO_error("write failed");
return -1;
}
);
default:
posthocIO_error("Unhandled data type");
return -1;
}
return 0;
}
}
namespace sensei
{
//-----------------------------------------------------------------------------
vtkStandardNewMacro(PosthocIO);
//-----------------------------------------------------------------------------
PosthocIO::PosthocIO() : Comm(MPI_COMM_WORLD), CommRank(0),
OutputDir("./"), HeaderFile("ImageHeader"), BlockExt(".sensei"),
HaveHeader(true), Mode(MpiIO), Period(1) {}
//-----------------------------------------------------------------------------
PosthocIO::~PosthocIO()
{}
//-----------------------------------------------------------------------------
void PosthocIO::Initialize(MPI_Comm comm,
const std::string &outputDir, const std::string &headerFile,
const std::string &blockExt, const std::vector<std::string> &cellArrays,
const std::vector<std::string> &pointArrays, int mode, int period)
{
#ifdef PosthocIO_DEBUG
posthocIO_status((this->CommRank==0), "PosthocIO::Initialize");
#endif
this->Comm = comm;
this->CommRank = 0;
MPI_Comm_rank(this->Comm, &this->CommRank);
this->OutputDir = outputDir;
this->HeaderFile = headerFile;
this->BlockExt = blockExt;
this->CellArrays = cellArrays;
this->PointArrays = pointArrays;
this->HaveHeader = (this->CommRank==0 ? false : true);
this->Mode = mode;
this->Period = period;
}
//-----------------------------------------------------------------------------
int PosthocIO::WriteBOVHeader(const std::string &fileName,
const std::vector<std::string> &arrays, const int *wholeExtent)
{
std::ofstream ff(fileName, std::ofstream::out);
if (!ff.good())
{
posthocIO_error("Failed to write the header file \"" << fileName << "\"")
return -1;
}
int dims[3] = {wholeExtent[1] - wholeExtent[0] + 1,
wholeExtent[3] - wholeExtent[2] + 1,
wholeExtent[5] - wholeExtent[4] + 1};
ff << "# SciberQuest MPI-IO BOV Reader" << std::endl
<< "nx=" << dims[0] << ", ny=" << dims[1] << ", nz=" << dims[2] << std::endl
<< "ext=" << this->BlockExt << std::endl
<< "dtype=f32" << std::endl
<< std::endl;
size_t n = arrays.size();
for (size_t i = 0; i < n; ++i)
ff << "scalar:" << arrays[i] << std::endl;
ff << std::endl;
ff.close();
#ifdef posthocIO_DEBUG
posthocIO_status((this->CommRank==0),
"wrote BOV header \"" << fileName << "\"");
#endif
return 0;
}
//-----------------------------------------------------------------------------
int PosthocIO::WriteBOVHeader(vtkInformation *info)
{
if (this->CommRank || this->HaveHeader)
return 0;
// handle both cell and point data
for (int dType = 0; dType < 2; ++dType)
{
// get the arrays
std::vector<std::string> &arrays =
dType ? this->CellArrays : this->PointArrays;
if (arrays.empty())
continue;
// get the extents
int wholeExt[6];
if (dType)
::getWholeCellExtents(info, wholeExt);
else
::getWholePointExtents(info, wholeExt);
const char *dTypeId =
(dType ? "CellData.bov" : "PointData.bov");
std::string headerFile =
this->OutputDir + "/" + this->HeaderFile + dTypeId;
this->WriteBOVHeader(headerFile, arrays, wholeExt);
}
this->HaveHeader = true;
return 0;
}
//-----------------------------------------------------------------------------
int PosthocIO::WriteXMLP(vtkCompositeDataSet *cd,
vtkInformation *info, int timeStep)
{
#ifdef PosthocIO_DEBUG
posthocIO_status((this->CommRank==0), "PosthocIO::WriteXMLP");
#endif
return 0;
}
//-----------------------------------------------------------------------------
int PosthocIO::WriteBOV(vtkCompositeDataSet *cd,
vtkInformation *info, int timeStep)
{
#ifdef PosthocIO_DEBUG
posthocIO_status((this->CommRank==0), "PosthocIO::WriteBOV");
#endif
// handle both cell and point data
for (int dType = 0; dType < 2; ++dType)
{
std::vector<std::string> &arrays =
dType ? this->CellArrays : this->PointArrays;
size_t n_arrays = arrays.size();
for (size_t i = 0; i < n_arrays; ++i)
{
const std::string &arrayName = arrays[i];
// construct the array's file name
std::ostringstream oss;