Updates will be applied April 15th at 12pm EDT (UTC-0400). GitLab could be a little slow between 12 - 12:45pm EDT.

Commit f56fe0f9 authored by Brad Whitlock's avatar Brad Whitlock Committed by Burlen Loring

Ghost cell/node support

Add capability to detect and add ghost masking arrays to the data
adaptor. The mask arrays follow VisIt's bit convention which is
now recognized by both ParaView and VisIt and use VTK's naming
convention. Full details can be found in the following links
https://blog.kitware.com/ghost-and-blanking-visibility-changes/
http://www.visitusers.org/index.php?title=Representing_ghost_data
The oscilators mini-apps has been updated as have the Catalist,
ADIOS, Libsim, autocorrelation and histogram analyses.
Co-Authored-By: Andrew Bauer's avatarAndy Bauer <andy.bauer@kitware.com>
Co-Authored-By: Burlen Loring's avatarBurlen Loring <bloring@lbl.gov>
parent 5ff37a10
......@@ -2,13 +2,17 @@
set enabled="1" on analyses you wish to enable -->
<sensei>
<!-- Custom Analyses-->
<analysis type="PosthocIO" array="data" association="cell"
enabled="0" period="1" mode="vtkXmlP" output_dir="./" />
<analysis type="PosthocIO"
output_dir="./" file_name="output" mode="visit" enabled="1">
<mesh name="mesh">
<cell_arrays> data </cell_arrays>
</mesh>
</analysis>
<analysis type="histogram" array="data" association="cell"
<analysis type="histogram" mesh="mesh" array="data" association="cell"
bins="10" enabled="1" />
<analysis type="autocorrelation" array="data" association="cell" window="10"
<analysis type="autocorrelation" mesh="mesh" array="data" association="cell" window="10"
k-max="3" enabled="1" />
<!-- VTK-m Analyses -->
......
......@@ -23,12 +23,13 @@ void initialize(MPI_Comm world,
int* gid,
int* from_x, int* from_y, int* from_z,
int* to_x, int* to_y, int* to_z,
int* shape, int ghostLevels,
const std::string& config_file)
{
(void)window;
timer::MarkEvent mark("oscillators::bridge::initialize");
GlobalDataAdaptor = vtkSmartPointer<oscillators::DataAdaptor>::New();
GlobalDataAdaptor->Initialize(nblocks);
GlobalDataAdaptor->Initialize(nblocks, shape, ghostLevels);
GlobalDataAdaptor->SetDataTimeStep(-1);
for (size_t cc=0; cc < n_local_blocks; ++cc)
......
......@@ -14,6 +14,7 @@ namespace bridge
int* gid,
int* from_x, int* from_y, int* from_z,
int* to_x, int* to_y, int* to_z,
int* shape, int ghostLevels,
const std::string& config_file);
void set_data(int gid, float* data);
......
......@@ -5,6 +5,7 @@
#include <vtkCellData.h>
#include <vtkFloatArray.h>
#include <vtkIdTypeArray.h>
#include <vtkIntArray.h>
#include <vtkImageData.h>
#include <vtkInformation.h>
#include <vtkMultiBlockDataSet.h>
......@@ -28,6 +29,8 @@ struct DataAdaptor::DInternals
vtkSmartPointer<vtkMultiBlockDataSet> uMesh;
std::vector<vtkSmartPointer<vtkUnstructuredGrid> > UnstructuredMesh;
std::vector<int> DataExtent;
int shape[3];
int ghostLevels;
};
inline bool areBoundsValid(const diy::DiscreteBounds& bds)
......@@ -52,7 +55,7 @@ DataAdaptor::~DataAdaptor()
}
//-----------------------------------------------------------------------------
void DataAdaptor::Initialize(size_t nblocks)
void DataAdaptor::Initialize(size_t nblocks, const int *shape_, int ghostLevels_)
{
DInternals& internals = (*this->Internals);
internals.CellExtents.resize(nblocks);
......@@ -68,6 +71,12 @@ void DataAdaptor::Initialize(size_t nblocks)
internals.CellExtents[cc].max[1] = -1;
internals.CellExtents[cc].max[2] = -1;
}
internals.shape[0] = shape_[0];
internals.shape[1] = shape_[1];
internals.shape[2] = shape_[2];
internals.ghostLevels = ghostLevels_;
this->ReleaseData();
}
......@@ -122,6 +131,8 @@ int DataAdaptor::GetMesh(const std::string &meshName, bool structureOnly,
{
internals.uMesh = vtkSmartPointer<vtkMultiBlockDataSet>::New();
internals.uMesh->SetNumberOfBlocks(static_cast<unsigned int>(internals.CellExtents.size()));
for (size_t cc=0; cc < internals.CellExtents.size(); ++cc)
internals.uMesh->SetBlock(static_cast<unsigned int>(cc), NULL);
}
// Either create empty vtkUnstructuredGrid objects or let us replace
// empty ones with new ones that have the right data.
......@@ -130,12 +141,20 @@ int DataAdaptor::GetMesh(const std::string &meshName, bool structureOnly,
unsigned int bid = static_cast<unsigned int>(cc);
vtkUnstructuredGrid *g = vtkUnstructuredGrid::SafeDownCast(
internals.uMesh->GetBlock(bid));
if((structureOnly && g == NULL) ||
(!structureOnly && g != NULL && g->GetNumberOfCells() == 0))
if(g == NULL)
{
g = (vtkUnstructuredGrid *)this->GetUnstructuredMesh(cc, structureOnly);
//cout << "Setting uMesh[" << bid << "] structureOnly=" << structureOnly << ", g=" << (void*)g << ", ncells=" << (g ? g->GetNumberOfCells() : 0) << endl;
internals.uMesh->SetBlock(bid, g);
}
else if(!structureOnly && g->GetNumberOfCells() == 0)
{
internals.uMesh->SetBlock(bid, this->GetUnstructuredMesh(cc, structureOnly));
g = (vtkUnstructuredGrid *)this->GetUnstructuredMesh(cc, structureOnly);
//cout << "Replacing uMesh[" << bid << "] structureOnly=" << structureOnly << ", g=" << (void*)g << ", ncells=" << (g ? g->GetNumberOfCells() : 0) << endl;
internals.uMesh->SetBlock(bid, g);
}
}
mesh = internals.uMesh;
}
else
......@@ -271,7 +290,7 @@ int DataAdaptor::AddArray(vtkDataObject* mesh, const std::string &meshName,
vtkMultiBlockDataSet* md = vtkMultiBlockDataSet::SafeDownCast(mesh);
for (unsigned int cc=0, max=md->GetNumberOfBlocks(); cc < max; ++cc)
{
if (!internals.Data[cc])
if (!internals.Data[cc]) // Exclude NULL datasets
{
continue;
}
......@@ -281,13 +300,13 @@ int DataAdaptor::AddArray(vtkDataObject* mesh, const std::string &meshName,
{
vtkSmartPointer<vtkImageData>& blockMesh = internals.BlockMesh[cc];
cd = (blockMesh? blockMesh->GetCellData() : NULL);
ncells = blockMesh->GetNumberOfCells();
ncells = (blockMesh? blockMesh->GetNumberOfCells() : 0);
}
else if(meshName == "ucdmesh")
{
vtkSmartPointer<vtkUnstructuredGrid>& uMesh = internals.UnstructuredMesh[cc];
cd = (uMesh? uMesh->GetCellData() : NULL);
ncells = uMesh->GetNumberOfCells();
ncells = (uMesh? uMesh->GetNumberOfCells() : 0);
}
if (cd != NULL)
{
......@@ -306,6 +325,142 @@ int DataAdaptor::AddArray(vtkDataObject* mesh, const std::string &meshName,
return retVal;
}
//----------------------------------------------------------------------------
int DataAdaptor::GetMeshHasGhostCells(const std::string &/*meshName*/,
int &nLayers)
{
DInternals& internals = (*this->Internals);
nLayers = internals.ghostLevels;
return 0;
}
//----------------------------------------------------------------------------
vtkDataArray *
DataAdaptor::CreateGhostCellsArray(int cc) const
{
// This sim is always 3D.
const DInternals& internals = (*this->Internals);
int imin = internals.CellExtents[cc].min[0];
int jmin = internals.CellExtents[cc].min[1];
int kmin = internals.CellExtents[cc].min[2];
int imax = internals.CellExtents[cc].max[0];
int jmax = internals.CellExtents[cc].max[1];
int kmax = internals.CellExtents[cc].max[2];
int nx = imax - imin + 1;
int ny = jmax - jmin + 1;
int nz = kmax - kmin + 1;
int nxny = nx*ny;
int ncells = nx * ny *nz;
int ng = internals.ghostLevels;
#define GCTYPE unsigned char
#define GCVTKARRAY vtkUnsignedCharArray
GCVTKARRAY *g = GCVTKARRAY::New();
g->SetNumberOfTuples(ncells);
memset(g->GetVoidPointer(0), 0, sizeof(GCTYPE) * ncells);
g->SetName("vtkGhostType");
GCTYPE *gptr = (GCTYPE *)g->GetVoidPointer(0);
GCTYPE ghost = 1;
if(imin > 0)
{
// Set the low I faces to ghosts.
for(int k = 0; k < nz; ++k)
for(int j = 0; j < ny; ++j)
for(int i = 0; i < ng; ++i)
gptr[k * nxny + j*nx + i] = ghost;
}
if(imax < internals.shape[0]-1)
{
// Set the high I faces to ghosts.
for(int k = 0; k < nz; ++k)
for(int j = 0; j < ny; ++j)
for(int i = nx-ng; i < nx; ++i)
gptr[k * nxny + j*nx + i] = ghost;
}
if(jmin > 0)
{
// Set the low J faces to ghosts.
for(int k = 0; k < nz; ++k)
for(int j = 0; j < ng; ++j)
for(int i = 0; i < nx; ++i)
gptr[k * nxny + j*nx + i] = ghost;
}
if(jmax < internals.shape[1]-1)
{
// Set the high J faces to ghosts.
for(int k = 0; k < nz; ++k)
for(int j = ny-ng; j < ny; ++j)
for(int i = 0; i < nx; ++i)
gptr[k * nxny + j*nx + i] = ghost;
}
if(kmin > 0)
{
// Set the low K faces to ghosts.
for(int k = 0; k < ng; ++k)
for(int j = 0; j < ny; ++j)
for(int i = 0; i < nx; ++i)
gptr[k * nxny + j*nx + i] = ghost;
}
if(kmax < internals.shape[2]-1)
{
// Set the high K faces to ghosts.
for(int k = nz-ng; k < nz; ++k)
for(int j = 0; j < ny; ++j)
for(int i = 0; i < nx; ++i)
gptr[k * nxny + j*nx + i] = ghost;
}
return g;
}
//----------------------------------------------------------------------------
int DataAdaptor::AddGhostCellsArray(vtkDataObject *mesh, const std::string &meshName)
{
#ifndef NDEBUG
if (meshName != "mesh" && meshName != "ucdmesh")
{
SENSEI_ERROR("the miniapp provides meshes \"mesh\" and \"ucdmesh\".")
return 1;
}
#else
(void)mesh;
(void)meshName;
#endif
int retVal = 1;
DInternals& internals = (*this->Internals);
vtkMultiBlockDataSet* md = vtkMultiBlockDataSet::SafeDownCast(mesh);
for (unsigned int cc=0, max=md->GetNumberOfBlocks(); cc < max; ++cc)
{
vtkCellData *cd = NULL;
if(meshName == "mesh")
{
vtkSmartPointer<vtkImageData>& blockMesh = internals.BlockMesh[cc];
cd = (blockMesh? blockMesh->GetCellData() : NULL);
}
else if(meshName == "ucdmesh")
{
vtkSmartPointer<vtkUnstructuredGrid>& uMesh = internals.UnstructuredMesh[cc];
cd = (uMesh? uMesh->GetCellData() : NULL);
}
// NOTE: we could end up with a NULL cd if we encounter an empty mesh slot
// since the vtkMultiBlockDataSet is mostly empty.
if (cd != NULL)
{
if (cd->GetArray("vtkGhostType") == NULL)
{
vtkDataArray *g = CreateGhostCellsArray(cc);
cd->AddArray(g);
g->Delete();
}
retVal = 0;
}
}
return retVal;
}
//-----------------------------------------------------------------------------
int DataAdaptor::GetNumberOfMeshes(unsigned int &numMeshes)
{
......@@ -335,6 +490,7 @@ int DataAdaptor::GetMeshName(unsigned int id, std::string &meshName)
int DataAdaptor::GetNumberOfArrays(const std::string &meshName, int association,
unsigned int &numberOfArrays)
{
numberOfArrays = 0;
if ((meshName == "mesh" || meshName == "ucdmesh") && (association == vtkDataObject::CELL))
{
numberOfArrays = 1;
......
......@@ -3,6 +3,8 @@
#include <sensei/DataAdaptor.h>
class vtkDataArray;
namespace oscillators
{
......@@ -16,7 +18,7 @@ public:
///
/// This initializes the data adaptor. This must be called once per simulation run.
/// @param nblocks is the total number of blocks in the simulation run.
void Initialize(size_t nblocks);
void Initialize(size_t nblocks, const int *shape, int ghostLevels);
/// @brief Set the extents for local blocks.
void SetBlockExtent(int gid,
......@@ -45,6 +47,10 @@ public:
int GetArrayName(const std::string &meshName, int association,
unsigned int index, std::string &arrayName) override;
int GetMeshHasGhostCells(const std::string &meshName, int &nLayers) override;
int AddGhostCellsArray(vtkDataObject* mesh, const std::string &meshName) override;
int ReleaseData() override;
protected:
......@@ -53,6 +59,7 @@ protected:
vtkDataObject* GetBlockMesh(int gid);
vtkDataObject* GetUnstructuredMesh(int gid, bool structureOnly);
vtkDataArray* CreateGhostCellsArray(int cc) const;
private:
DataAdaptor(const DataAdaptor&); // not implemented.
......
......@@ -96,6 +96,7 @@ int main(int argc, char** argv)
size_t window = 10;
size_t k_max = 3;
int threads = 1;
int ghostLevels = 0;
std::string config_file;
std::string out_prefix = "";
Options ops(argc, argv);
......@@ -112,6 +113,7 @@ int main(int argc, char** argv)
>> Option( "t-end", t_end, "end time")
>> Option('j', "jobs", threads, "number of threads to use")
>> Option('o', "output", out_prefix, "prefix to save output")
>> Option('g', "ghost levels", ghostLevels, "Number of ghost levels")
;
bool sync = ops >> Present("sync", "synchronize after each time step");
bool log = ops >> Present("log", "generate full time and memory usage log");
......@@ -176,6 +178,9 @@ int main(int argc, char** argv)
to_x, to_y, to_z;
diy::RegularDecomposer<Bounds>::BoolVector share_face, wrap;
diy::RegularDecomposer<Bounds>::CoordinateVector ghosts = {ghostLevels, ghostLevels, ghostLevels};
// decompose the domain
diy::decompose(3, world.rank(), domain, assigner,
[&](int gid, const Bounds&, const Bounds& bounds, const Bounds& domain, const Link& link)
......@@ -192,7 +197,8 @@ int main(int argc, char** argv)
to_x.push_back(bounds.max[0]);
to_y.push_back(bounds.max[1]);
to_z.push_back(bounds.max[2]);
});
},
share_face, wrap, ghosts);
timer::MarkEndEvent("oscillators::initialize");
......@@ -203,6 +209,7 @@ int main(int argc, char** argv)
&gids[0],
&from_x[0], &from_y[0], &from_z[0],
&to_x[0], &to_y[0], &to_z[0],
&shape[0], ghostLevels,
config_file);
#else
init_analysis(world, window, gids.size(),
......
......@@ -82,9 +82,10 @@ bool ADIOSAnalysisAdaptor::Execute(DataAdaptor* dataAdaptor)
SENSEI_WARNING("No subset specified. Writing all available data")
}
// collect the specified data objects
// collect the specified data objects and metadata
std::vector<vtkDataObject*> objects;
std::vector<std::string> objectNames;
std::vector<std::pair<int,int>> ghostLayers;
MeshRequirementsIterator mit =
this->Requirements.GetMeshRequirementsIterator();
......@@ -99,6 +100,42 @@ bool ADIOSAnalysisAdaptor::Execute(DataAdaptor* dataAdaptor)
return -1;
}
// get ghost cell/node metadata always provide this information as
// it is essential to process the data objects
int nGhostCellLayers = 0;
int nGhostNodeLayers = 0;
if (dataAdaptor->GetMeshHasGhostCells(mit.MeshName(), nGhostCellLayers) ||
dataAdaptor->GetMeshHasGhostNodes(mit.MeshName(), nGhostNodeLayers))
{
SENSEI_ERROR("Failed to get ghost layer info for mesh \"" << mit.MeshName() << "\"")
return -1;
}
// pass ghost layer metadata in field data.
vtkIntArray *glmd = vtkIntArray::New();
glmd->SetName("senseiGhostLayers");
glmd->SetNumberOfTuples(2);
glmd->SetValue(0, nGhostCellLayers);
glmd->SetValue(1, nGhostNodeLayers);
vtkFieldData *fd = dobj->GetFieldData();
fd->AddArray(glmd);
glmd->Delete();
// add the ghost cell arrays to the mesh
if ((nGhostCellLayers > 0) && dataAdaptor->AddGhostCellsArray(dobj, mit.MeshName()))
{
SENSEI_ERROR("Failed to get ghost cells for mesh \"" << mit.MeshName() << "\"")
return -1;
}
// add the ghost node arrays to the mesh
if (nGhostNodeLayers > 0 && dataAdaptor->AddGhostNodesArray(dobj, mit.MeshName()))
{
SENSEI_ERROR("Failed to get ghost nodes for mesh \"" << mit.MeshName() << "\"")
return -1;
}
// add the required arrays
ArrayRequirementsIterator ait =
this->Requirements.GetArrayRequirementsIterator(mit.MeshName());
......
......@@ -4,6 +4,7 @@
#include "Timer.h"
#include "ADIOSSchema.h"
#include "VTKUtils.h"
#include "MeshMetadata.h"
#include <vtkCompositeDataIterator.h>
#include <vtkDataSetAttributes.h>
......@@ -15,51 +16,46 @@
#include <sstream>
using vtkDataObjectPtr = vtkSmartPointer<vtkDataObject>;
// associates an object and flag indicating if it is static
// to a name
using ObjectMapType = std::map<std::string, std::pair<vtkDataObjectPtr, int>>;
namespace sensei
{
struct ADIOSDataAdaptor::InternalsType
// associate the mesh name to a data object and mesh metadata
using vtkDataObjectPtr = vtkSmartPointer<vtkDataObject>;
using ObjectType = std::pair<vtkDataObjectPtr, MeshMetadata>;
using ObjectMapType = std::map<std::string, ObjectType>;
using ObjectMapIterType = ObjectMapType::iterator;
// helpers to keep things legible accessing pairs of pairs
static void setObject(ObjectType &ot, vtkDataObject *dobj){ ot.first.TakeReference(dobj); }
static vtkDataObject *getObject(ObjectType &ot) { return ot.first.GetPointer(); }
static vtkDataObjectPtr &getObjectPtr(ObjectType &ot) { return ot.first; }
static MeshMetadata &getMetadata(ObjectType &ot){ return ot.second; }
static const std::string &getObjectName(ObjectMapIterType &it) { return it->first; }
static void setObject(ObjectMapIterType &it, vtkDataObject *dobj){ setObject(it->second, dobj); }
static vtkDataObject *getObject(ObjectMapIterType &it) { return getObject(it->second); }
static vtkDataObjectPtr &getObjectPtr(ObjectMapIterType &it) { return getObjectPtr(it->second); }
static MeshMetadata &getMetadata(ObjectMapIterType &it){ return getMetadata(it->second); }
static ObjectMapIterType find(ObjectMapType &objMap, const std::string &meshName)
{
InternalsType() : Comm(MPI_COMM_WORLD), Stream() {}
ObjectMapType::iterator FindObject(const std::string &name)
{
return this->ObjectMap.find(name);
}
vtkDataObject *GetObject(ObjectMapType::iterator &it)
{
return it->second.first.GetPointer();
}
return objMap.find(meshName);
}
const vtkDataObject *GetObject(ObjectMapType::const_iterator &it)
{
return it->second.first.GetPointer();
}
static bool good(ObjectMapType &objMap, const ObjectMapIterType &it)
{
return objMap.end() != it;
}
int GetStatic(ObjectMapType::iterator &it)
{
return it->second.second;
}
int GetStatic(ObjectMapType::const_iterator &it)
{
return it->second.second;
}
struct ADIOSDataAdaptor::InternalsType
{
InternalsType() : Comm(MPI_COMM_WORLD), Stream() {}
MPI_Comm Comm;
senseiADIOS::InputStream Stream;
senseiADIOS::DataObjectCollectionSchema Schema;
std::set<std::string> ObjectNames;
ObjectMapType ObjectMap;
};
......@@ -82,7 +78,7 @@ ADIOSDataAdaptor::~ADIOSDataAdaptor()
//----------------------------------------------------------------------------
void ADIOSDataAdaptor::EnableDynamicMesh(const std::string &meshName, int val)
{
this->Internals->ObjectMap[meshName].second = val;
getMetadata(this->Internals->ObjectMap[meshName]).StaticMesh = val;
}
//----------------------------------------------------------------------------
......@@ -140,7 +136,6 @@ int ADIOSDataAdaptor::Close()
timer::MarkEvent mark("ADIOSDataAdaptor::Close");
this->Internals->ObjectMap.clear();
this->Internals->ObjectNames.clear();
this->Internals->Stream.Close();
return 0;
......@@ -180,8 +175,6 @@ int ADIOSDataAdaptor::UpdateTimeStep()
this->SetDataTime(time);
// update the available meshes
this->Internals->ObjectNames.clear();
std::vector<std::string> names;
if (this->Internals->Schema.ReadObjectNames(this->Internals->Comm,
this->Internals->Stream, names))
......@@ -190,9 +183,14 @@ int ADIOSDataAdaptor::UpdateTimeStep()
return -1;
}
// try to add new object, if one exists nothing changes, preserving metadata
// this is necessary since static flag can be set before mesh names are read
unsigned int nNames = names.size();
for (unsigned int i = 0; i < nNames; ++i)
this->Internals->ObjectNames.insert(names[i]);
{
this->Internals->ObjectMap.insert(std::make_pair(names[i],
ObjectType(vtkDataObjectPtr(), MeshMetadata(names[i]))));
}
return 0;
}
......@@ -200,7 +198,7 @@ int ADIOSDataAdaptor::UpdateTimeStep()
//----------------------------------------------------------------------------
int ADIOSDataAdaptor::GetNumberOfMeshes(unsigned int &numMeshes)
{
numMeshes = this->Internals->ObjectNames.size();
numMeshes = this->Internals->ObjectMap.size();
return 0;
}
......@@ -209,121 +207,174 @@ int ADIOSDataAdaptor::GetMeshName(unsigned int id, std::string &meshName)
{
meshName = "";
if (id >= this->Internals->ObjectNames.size())
if (id >= this->Internals->ObjectMap.size())
{
SENSEI_ERROR("Mesh name " << id << " out of bounds. "
<< " only " << this->Internals->ObjectNames.size()
<< " only " << this->Internals->ObjectMap.size()
<< " mesh names available")
return -1;
}
std::set<std::string>::iterator it = this->Internals->ObjectNames.begin();
ObjectMapIterType it = this->Internals->ObjectMap.begin();
for (unsigned int i = 0; i < id; ++i)
++it;
meshName = *it;
meshName = getObjectName(it);
return 0;
}
//----------------------------------------------------------------------------
int ADIOSDataAdaptor::GetMesh(const std::string &meshName,
bool structure_only, vtkDataObject *&mesh)
bool structureOnly, vtkDataObject *&mesh)
{
timer::MarkEvent mark("ADIOSDataAdaptor::GetMesh");
mesh = nullptr;
int staticMesh = 1;
ObjectMapType::iterator it = this->Internals->FindObject(meshName);