Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-10-05 01:49:02 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Start creating test for HierarchicalHyperSweeper and fix various bugs in it

Browse Source
This commit is contained in:
Gunther H. Weber 2021-07-30 19:23:23 -07:00
parent 45a3570270
commit 3cf68a512e
21 changed files with 665 additions and 324 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
data/data/misc/8x9test_HierarchicalAugmentedTree_Block0.dat Normal file
View File

@ -0,0 +1,3 @@
version https://git-lfs.github.com/spec/v1
oid sha256:c5c3afc7bdb0fa75bad78fb69ce0557c752b98cab6f000d757ef2706a4498576
size 2116

3
data/data/misc/8x9test_HierarchicalAugmentedTree_Block1.dat Normal file
View File

@ -0,0 +1,3 @@
version https://git-lfs.github.com/spec/v1
oid sha256:e31bdb91419a32e98dc1061ae6ce348f8fbc6c79b26a4e0bcd352d2b2e62b6a8
size 2228

3
data/data/misc/8x9test_HierarchicalAugmentedTree_Block2.dat Normal file
View File

@ -0,0 +1,3 @@
version https://git-lfs.github.com/spec/v1
oid sha256:ccbe5767ce622b4c0e5e93dd6f6243d5c76a5a6f07c1e5cf43fb2b5e2e069255
size 2076

3
data/data/misc/8x9test_HierarchicalAugmentedTree_Block3.dat Normal file
View File

@ -0,0 +1,3 @@
version https://git-lfs.github.com/spec/v1
oid sha256:6e9f0913c17a3a12338d043c455c6e9cde51c72fa70d619d935fa4369effed45
size 2316

4
vtkm/worklet/contourtree_augmented/DataSetMesh.h
View File

@ -203,7 +203,7 @@ protected:
template <typename MeshTypeObj>
void GetOwnedVerticesByGlobalIdImpl(
const MeshTypeObj* mesh,
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler* localToGlobalIdRelabeler,
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler& localToGlobalIdRelabeler,
IdArrayType& ownedVertices) const;
virtual void DebugPrintExtends();
@ -215,7 +215,7 @@ protected:
template <typename MeshTypeObj>
void DataSetMesh::GetOwnedVerticesByGlobalIdImpl(
const MeshTypeObj* mesh,
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler* localToGlobalIdRelabeler,
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler& localToGlobalIdRelabeler,
IdArrayType& ownedVertices) const
{
// use temporary array since we need to compress these at the end via CopyIf so we

4
vtkm/worklet/contourtree_augmented/data_set_mesh/GetOwnedVerticesByGlobalIdWorklet.h
View File

@ -78,7 +78,7 @@ public:
// Constructor
VTKM_EXEC_CONT
GetOwnedVerticesByGlobalIdWorklet(
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler* localToGlobalIdRelabeler)
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler& localToGlobalIdRelabeler)
: LocalToGlobalIdRelabeler(localToGlobalIdRelabeler)
{
}
@ -93,7 +93,7 @@ public:
}
private:
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler* LocalToGlobalIdRelabeler;
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler& LocalToGlobalIdRelabeler;
}; // Mesh2D_DEM_VertexStarter
} // namespace data_set_mesh

12
vtkm/worklet/contourtree_augmented/data_set_mesh/MeshStructure2D.h
View File

@ -101,7 +101,7 @@ public:
/// otherwise it returns global id of the vertex as determined via the IdRelabeler
VTKM_EXEC_CONT
inline vtkm::Id GetVertexOwned(const vtkm::Id& meshIndex,
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler*
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler&
localToGlobalIdRelabeler) const
{
// Get the vertex position
@ -109,17 +109,17 @@ public:
// now test - the low ID boundary belongs to this block
// the high ID boundary belongs to the next block if there is one
if (((pos[1] == this->MeshSize[1] - 1) &&
(pos[1] + localToGlobalIdRelabeler->LocalBlockOrigin[1] !=
localToGlobalIdRelabeler->GlobalSize[1] - 1)) ||
(pos[1] + localToGlobalIdRelabeler.LocalBlockOrigin[1] !=
localToGlobalIdRelabeler.GlobalSize[1] - 1)) ||
((pos[0] == this->MeshSize[0] - 1) &&
(pos[0] + localToGlobalIdRelabeler->LocalBlockOrigin[0] !=
localToGlobalIdRelabeler->GlobalSize[0] - 1)))
(pos[0] + localToGlobalIdRelabeler.LocalBlockOrigin[0] !=
localToGlobalIdRelabeler.GlobalSize[0] - 1)))
{
return vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT;
}
else
{
return (*localToGlobalIdRelabeler)(meshIndex);
return localToGlobalIdRelabeler(meshIndex);
}
}

17
vtkm/worklet/contourtree_augmented/data_set_mesh/MeshStructure3D.h
View File

@ -109,7 +109,7 @@ public:
/// otherwise it returns global id of the vertex as determined via the IdRelabeler
VTKM_EXEC_CONT
inline vtkm::Id GetVertexOwned(const vtkm::Id& meshIndex,
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler*
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler&
localToGlobalIdRelabeler) const
{
// Get the vertex position
@ -117,24 +117,23 @@ public:
// now test - the low ID boundary belongs to this block
// the high ID boundary belongs to the next block if there is one
if (((pos[1] == this->MeshSize[1] - 1) &&
(pos[1] + localToGlobalIdRelabeler->LocalBlockOrigin[1] !=
localToGlobalIdRelabeler->GlobalSize[1] - 1)) ||
(pos[1] + localToGlobalIdRelabeler.LocalBlockOrigin[1] !=
localToGlobalIdRelabeler.GlobalSize[1] - 1)) ||
((pos[0] == this->MeshSize[0] - 1) &&
(pos[0] + localToGlobalIdRelabeler->LocalBlockOrigin[0] !=
localToGlobalIdRelabeler->GlobalSize[0] - 1)) ||
(pos[0] + localToGlobalIdRelabeler.LocalBlockOrigin[0] !=
localToGlobalIdRelabeler.GlobalSize[0] - 1)) ||
((pos[2] == this->MeshSize[2] - 1) &&
(pos[2] + localToGlobalIdRelabeler->LocalBlockOrigin[2] !=
localToGlobalIdRelabeler->GlobalSize[2] - 1)))
(pos[2] + localToGlobalIdRelabeler.LocalBlockOrigin[2] !=
localToGlobalIdRelabeler.GlobalSize[2] - 1)))
{
return vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT;
}
else
{
return (*localToGlobalIdRelabeler)(meshIndex);
return localToGlobalIdRelabeler(meshIndex);
}
}
vtkm::Id3 MeshSize;
}; // Mesh_DEM_2D_ExecutionObject

2
vtkm/worklet/contourtree_augmented/data_set_mesh/NotNoSuchElementPredicate.h
View File

@ -75,7 +75,7 @@ public:
VTKM_EXEC_CONT
bool operator()(const vtkm::Id& meshVertexId) const
{
return static_cast<bool>(vtkm::worklet::contourtree_augmented::NoSuchElement(meshVertexId));
return bool{ !vtkm::worklet::contourtree_augmented::NoSuchElement(meshVertexId) };
}
private:

6
vtkm/worklet/contourtree_augmented/meshtypes/DataSetMeshTriangulation2DFreudenthal.h
View File

@ -109,7 +109,7 @@ public:
/// Get of global indices of the vertices owned by this mesh. Implemented via
/// DataSetMesh.GetOwnedVerticesByGlobalIdImpl
void GetOwnedVerticesByGlobalId(
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler* localToGlobalIdRelabeler,
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler& localToGlobalIdRelabeler,
IdArrayType& ownedVertices) const;
private:
@ -175,8 +175,8 @@ inline void DataSetMeshTriangulation2DFreudenthal::GetBoundaryVertices(
}
// Overwrite the implemenation from the base DataSetMesh parent class
void DataSetMeshTriangulation2DFreudenthal::GetOwnedVerticesByGlobalId(
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler* localToGlobalIdRelabeler,
inline void DataSetMeshTriangulation2DFreudenthal::GetOwnedVerticesByGlobalId(
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler& localToGlobalIdRelabeler,
IdArrayType& ownedVertices) const
{
return this->GetOwnedVerticesByGlobalIdImpl(this, localToGlobalIdRelabeler, ownedVertices);

6
vtkm/worklet/contourtree_augmented/meshtypes/DataSetMeshTriangulation3DFreudenthal.h
View File

@ -103,7 +103,7 @@ public:
/// Get of global indices of the vertices owned by this mesh. Implemented via
/// DataSetMesh.GetOwnedVerticesByGlobalIdImpl.
void GetOwnedVerticesByGlobalId(
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler* localToGlobalIdRelabeler,
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler& localToGlobalIdRelabeler,
IdArrayType& ownedVertices) const;
private:
@ -180,8 +180,8 @@ inline void DataSetMeshTriangulation3DFreudenthal::GetBoundaryVertices(
}
// Overwrite the implemenation from the base DataSetMesh parent class
void DataSetMeshTriangulation3DFreudenthal::GetOwnedVerticesByGlobalId(
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler* localToGlobalIdRelabeler,
inline void DataSetMeshTriangulation3DFreudenthal::GetOwnedVerticesByGlobalId(
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler& localToGlobalIdRelabeler,
IdArrayType& ownedVertices) const
{
return this->GetOwnedVerticesByGlobalIdImpl(this, localToGlobalIdRelabeler, ownedVertices);

6
vtkm/worklet/contourtree_augmented/meshtypes/DataSetMeshTriangulation3DMarchingCubes.h
View File

@ -108,7 +108,7 @@ public:
/// Get of global indices of the vertices owned by this mesh. Implemented via
/// DataSetMesh.GetOwnedVerticesByGlobalIdImpl.
void GetOwnedVerticesByGlobalId(
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler* localToGlobalIdRelabeler,
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler& localToGlobalIdRelabeler,
IdArrayType& ownedVertices) const;
private:
@ -213,8 +213,8 @@ inline void DataSetMeshTriangulation3DMarchingCubes::GetBoundaryVertices(
}
// Overwrite the implemenation from the base DataSetMesh parent class
void DataSetMeshTriangulation3DMarchingCubes::GetOwnedVerticesByGlobalId(
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler* localToGlobalIdRelabeler,
inline void DataSetMeshTriangulation3DMarchingCubes::GetOwnedVerticesByGlobalId(
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler& localToGlobalIdRelabeler,
IdArrayType& ownedVertices) const
{
return this->GetOwnedVerticesByGlobalIdImpl(this, localToGlobalIdRelabeler, ownedVertices);

61
vtkm/worklet/contourtree_distributed/HierarchicalContourTree.h
View File

@ -71,8 +71,10 @@
#define VOLUME_PRINT_WIDTH 8
#include <vtkm/Types.h>
#include <vtkm/worklet/contourtree_augmented/ContourTree.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/ContourTreeMesh.h>
#include <vtkm/worklet/contourtree_distributed/LoadArrays.h>
#include <vtkm/worklet/contourtree_distributed/hierarchical_contour_tree/FindRegularByGlobal.h>
#include <vtkm/worklet/contourtree_distributed/hierarchical_contour_tree/FindSuperArcForUnknownNode.h>
#include <vtkm/worklet/contourtree_distributed/hierarchical_contour_tree/InitalizeSuperchildrenWorklet.h>
@ -143,8 +145,8 @@ public:
// how many rounds of fan-in were used to construct it
vtkm::Id NumRounds;
// use for debugging?
vtkm::Id NumOwnedRegularVertices;
// use for debugging? -> This makes more sense in hyper sweeper?
// vtkm::Id NumOwnedRegularVertices;
// The following arrays store the numbers of reg/super/hyper nodes at each level of the hierarchy
// They are filled in from the top down, and are fundamentally CPU side control variables
@ -165,7 +167,7 @@ public:
/// routine to create a FindRegularByGlobal object that we can use as an input for worklets to call the function
VTKM_CONT
FindRegularByGlobal GetFindRegularByGlobal()
FindRegularByGlobal GetFindRegularByGlobal() const
{
return FindRegularByGlobal(this->RegularNodeSortOrder, this->RegularNodeGlobalIds);
}
@ -215,19 +217,24 @@ public:
VTKM_CONT
std::string PrintDotSuperStructure(const char* label) const;
/// Load from file saved by PPP
VTKM_CONT
void Load(const char* filename);
/// Print hierarchical tree construction stats, usually used for logging
VTKM_CONT
std::string PrintTreeStats() const;
/// debug routine
VTKM_CONT
std::string DebugPrint(const char* message, const char* fileName, long lineNum) const;
std::string DebugPrint(std::string message, const char* fileName, long lineNum) const;
// modified version of dumpSuper() that also gives volume counts
VTKM_CONT
std::string DumpVolumes(vtkm::Id totalVolume,
const vtkm::worklet::contourtree_augmented::IdArrayType& intrinsicVolume,
const vtkm::worklet::contourtree_augmented::IdArrayType& dependentVolume);
std::string DumpVolumes(
vtkm::Id totalVolume,
const vtkm::worklet::contourtree_augmented::IdArrayType& intrinsicVolume,
const vtkm::worklet::contourtree_augmented::IdArrayType& dependentVolume) const;
private:
/// Used internally to Invoke worklets
@ -236,7 +243,7 @@ private:
template <typename FieldType>
HierarchicalContourTree<FieldType>::HierarchicalContourTree()
: NumOwnedRegularVertices(static_cast<vtkm::Id>(0))
//: NumOwnedRegularVertices(static_cast<vtkm::Id>(0))
{ // constructor
NumRegularNodesInRound.ReleaseResources();
NumSupernodesInRound.ReleaseResources();
@ -749,9 +756,41 @@ std::string HierarchicalContourTree<FieldType>::PrintDotSuperStructure(const cha
return std::string("HierarchicalContourTree<FieldType>::PrintDotSuperStructure() Complete");
} // PrintDotSuperStructure
template <typename FieldType>
void HierarchicalContourTree<FieldType>::Load(const char* filename)
{
std::ifstream is(filename);
ReadIndexArray(is, this->RegularNodeGlobalIds);
ReadDataArray<FieldType>(is, this->DataValues);
ReadIndexArray(is, this->RegularNodeSortOrder);
ReadIndexArray(is, this->Regular2Supernode);
ReadIndexArray(is, this->Superparents);
ReadIndexArray(is, this->Supernodes);
ReadIndexArray(is, this->Superarcs);
ReadIndexArray(is, this->Hyperparents);
ReadIndexArray(is, this->Super2Hypernode);
ReadIndexArray(is, this->WhichRound);
ReadIndexArray(is, this->WhichIteration);
ReadIndexArray(is, this->Hypernodes);
ReadIndexArray(is, this->Hyperarcs);
ReadIndexArray(is, this->Superchildren);
int nRounds;
is.read(reinterpret_cast<char*>(&nRounds), sizeof(nRounds));
//std::cout << "nRounds = " << nRounds << std::endl;
this->NumRounds = nRounds;
//this->NumOwnedRegularVertices = 0;
ReadIndexArray(is, this->NumRegularNodesInRound);
ReadIndexArray(is, this->NumSupernodesInRound);
ReadIndexArray(is, this->NumHypernodesInRound);
ReadIndexArray(is, this->NumIterations);
ReadIndexArrayVector(is, this->FirstSupernodePerIteration);
ReadIndexArrayVector(is, this->FirstHypernodePerIteration);
}
/// debug routine
template <typename FieldType>
std::string HierarchicalContourTree<FieldType>::DebugPrint(const char* message,
std::string HierarchicalContourTree<FieldType>::DebugPrint(std::string message,
const char* fileName,
long lineNum) const
{ // DebugPrint
@ -811,7 +850,7 @@ std::string HierarchicalContourTree<FieldType>::DebugPrint(const char* message,
"nSupernodes In Round", this->NumSupernodesInRound, -1, resultStream);
vtkm::worklet::contourtree_augmented::PrintIndices(
"nHypernodes In Round", this->NumHypernodesInRound, -1, resultStream);
resultStream << "Owned Regular Vertices: " << this->NumOwnedRegularVertices << std::endl;
//resultStream << "Owned Regular Vertices: " << this->NumOwnedRegularVertices << std::endl;
vtkm::worklet::contourtree_augmented::PrintHeader(this->NumIterations.GetNumberOfValues(),
resultStream);
vtkm::worklet::contourtree_augmented::PrintIndices(
@ -861,7 +900,7 @@ template <typename FieldType>
std::string HierarchicalContourTree<FieldType>::DumpVolumes(
vtkm::Id totalVolume,
const vtkm::worklet::contourtree_augmented::IdArrayType& intrinsicVolume,
const vtkm::worklet::contourtree_augmented::IdArrayType& dependentVolume)
const vtkm::worklet::contourtree_augmented::IdArrayType& dependentVolume) const
{ // DumpVolumes()
// a local string stream to build the output
std::stringstream outStream;

354
vtkm/worklet/contourtree_distributed/HierarchicalHyperSweeper.h
View File

@ -102,15 +102,12 @@ namespace contourtree_distributed
{
// the class itself
template <typename MeshType, typename FieldType>
template <typename SweepValueType, typename ContourTreeFieldType>
class HierarchicalHyperSweeper
{ // class HierarchicalHyperSweeper
public:
// the tree that it hypersweeps over
const HierarchicalContourTree<FieldType>& HierarchicalTree;
// the underlying mesh base block type
const MeshType& BaseBlock;
const HierarchicalContourTree<ContourTreeFieldType>& HierarchicalTree;
// the Id of the base block (used for debug output)
vtkm::Id BlockId;
@ -118,9 +115,9 @@ public:
// array of values being operated over (same size as supernode set)
// keep both intrinsic & dependent values
// the intrinsic values are just stored but not modifid here
const vtkm::cont::ArrayHandle<FieldType>& IntrinsicValues;
const vtkm::cont::ArrayHandle<SweepValueType>& IntrinsicValues;
// the dependent values are what is being sweeped and are updated here
const vtkm::cont::ArrayHandle<FieldType>& DependentValues;
const vtkm::cont::ArrayHandle<SweepValueType>& DependentValues;
// and to avoid an extra log summation, store the number of logical nodes for the underlying block
// (computed when initializing the regular vertex list)
vtkm::Id NumOwnedRegularVertices;
@ -129,7 +126,7 @@ public:
// these are working arrays, lifted up here for ease of debug code
// Subranges of these arrays will be reused in the rounds / iterations rather than being reallocated
// an array for temporary storage of the prefix sums
vtkm::cont::ArrayHandle<FieldType> ValuePrefixSum;
vtkm::cont::ArrayHandle<SweepValueType> ValuePrefixSum;
// two arrays for collecting targets of transfers
vtkm::worklet::contourtree_augmented::IdArrayType TransferTarget;
vtkm::worklet::contourtree_augmented::IdArrayType SortedTransferTarget;
@ -142,12 +139,11 @@ public:
/// @param[in] baseBlock the underlying mesh base block type
/// @param[in] intrinsicValues array of values of intrinisic nodes are just being stored here but not modified
/// @param[in] dependentValues array of values being operated over (same size as supernode set)
HierarchicalHyperSweeper<MeshType, FieldType>(
HierarchicalHyperSweeper<SweepValueType, ContourTreeFieldType>(
vtkm::Id blockId,
const HierarchicalContourTree<FieldType>& hierarchicalTree,
const MeshType& baseBlock,
const vtkm::cont::ArrayHandle<FieldType>& intrinsicValues,
const vtkm::cont::ArrayHandle<FieldType>& dependentValues);
const HierarchicalContourTree<ContourTreeFieldType>& hierarchicalTree,
const vtkm::cont::ArrayHandle<SweepValueType>& intrinsicValues,
const vtkm::cont::ArrayHandle<SweepValueType>& dependentValues);
/// Our routines to initialize the sweep need to be static (or externa)l if we are going to use the constructor
/// to run the actual hypersweep
@ -155,10 +151,11 @@ public:
/// @param[in] baseBlock the underlying mesh base block to initialize from
/// @param[in] localToGlobalIdRelabeler Id relabeler used to compute global indices from local mesh indices
/// @param[out] superarcRegularCounts arrray for the output superarc regular counts
static void InitializeIntrinsicVertexCount(
const HierarchicalContourTree<FieldType>& hierarchicalTree,
template <typename MeshType>
void InitializeIntrinsicVertexCount(
const HierarchicalContourTree<ContourTreeFieldType>& hierarchicalTree,
const MeshType& baseBlock,
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler* localToGlobalIdRelabeler,
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler& localToGlobalIdRelabeler,
vtkm::worklet::contourtree_augmented::IdArrayType& superarcRegularCounts);
/// routine to do the local hypersweep using addition / subtraction
@ -205,27 +202,27 @@ private:
}; // class HierarchicalHyperSweeper
template <typename MeshType, typename FieldType>
HierarchicalHyperSweeper<MeshType, FieldType>::HierarchicalHyperSweeper(
template <typename SweepValueType, typename ContourTreeFieldType>
HierarchicalHyperSweeper<SweepValueType, ContourTreeFieldType>::HierarchicalHyperSweeper(
vtkm::Id blockId,
const HierarchicalContourTree<FieldType>& hierarchicalTree,
const MeshType& baseBlock,
const vtkm::cont::ArrayHandle<FieldType>& intrinsicValues,
const vtkm::cont::ArrayHandle<FieldType>& dependentValues)
const HierarchicalContourTree<ContourTreeFieldType>& hierarchicalTree,
const vtkm::cont::ArrayHandle<SweepValueType>& intrinsicValues,
const vtkm::cont::ArrayHandle<SweepValueType>& dependentValues)
: HierarchicalTree(hierarchicalTree)
, BaseBlock(baseBlock)
, BlockId(blockId)
, IntrinsicValues(intrinsicValues)
, DependentValues(dependentValues)
, NumOwnedRegularVertices(static_cast<vtkm::Id>(0))
, NumOwnedRegularVertices(vtkm::Id{ 0 })
{ // constructor
// Initalize arrays with 0s
vtkm::cont::ArrayHandleConstant<vtkm::Id> tempZeroArray(
0, this->HierarchicalTree.Supernodes.GetNumberOfValues());
vtkm::cont::Algorithm::Copy(tempZeroArray, this->ValuePrefixSum);
vtkm::cont::Algorithm::Copy(tempZeroArray, this->TransferTarget);
vtkm::cont::Algorithm::Copy(tempZeroArray, this->SortedTransferTarget);
// initialise the supersortPermute to the identity
vtkm::cont::Algorithm::Fill(
this->ValuePrefixSum, vtkm::Id{ 0 }, this->HierarchicalTree.Supernodes.GetNumberOfValues());
vtkm::cont::Algorithm::Fill(
this->TransferTarget, vtkm::Id{ 0 }, this->HierarchicalTree.Supernodes.GetNumberOfValues());
vtkm::cont::Algorithm::Fill(this->SortedTransferTarget,
vtkm::Id{ 0 },
this->HierarchicalTree.Supernodes.GetNumberOfValues());
// Initialize the supersortPermute to the identity
vtkm::cont::ArrayHandleIndex tempIndexArray(
this->HierarchicalTree.Supernodes.GetNumberOfValues());
vtkm::cont::Algorithm::Copy(tempIndexArray, this->SuperSortPermute);
@ -233,15 +230,14 @@ HierarchicalHyperSweeper<MeshType, FieldType>::HierarchicalHyperSweeper(
// static function used to compute the initial superarc regular counts
template <typename MeshType, typename FieldType>
void HierarchicalHyperSweeper<MeshType, FieldType>::InitializeIntrinsicVertexCount(
const HierarchicalContourTree<FieldType>& hierarchicalTree,
template <typename SweepValueType, typename ContourTreeFieldType>
template <typename MeshType>
void HierarchicalHyperSweeper<SweepValueType, ContourTreeFieldType>::InitializeIntrinsicVertexCount(
const HierarchicalContourTree<ContourTreeFieldType>& hierarchicalTree,
const MeshType& baseBlock,
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler* localToGlobalIdRelabeler,
const vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler& localToGlobalIdRelabeler,
vtkm::worklet::contourtree_augmented::IdArrayType& superarcRegularCounts)
{ // InitializeIntrinsicVertexCount()
vtkm::cont::Invoker
localInvoke; // Needed because this a static function so we can't use the invoke from the object
// I. Call the mesh to get a list of all regular vertices belonging to the block by global Id
vtkm::worklet::contourtree_augmented::IdArrayType globalIds;
// TODO/FIXME: Even though the virtual function on DataSetMesh was removed in commit
@ -250,14 +246,15 @@ void HierarchicalHyperSweeper<MeshType, FieldType>::InitializeIntrinsicVertexCou
// this is indeed correct.
baseBlock.GetOwnedVerticesByGlobalId(localToGlobalIdRelabeler, globalIds);
// and store the size for later reference
hierarchicalTree.NumOwnedRegularVertices = globalIds.GetNumberOfValues();
//hierarchicalTree.NumOwnedRegularVertices = globalIds.GetNumberOfValues();
this->NumOwnedRegularVertices = globalIds.GetNumberOfValues();
#ifdef DEBUG_PRINT
{
std::stringstream debugStream;
debugStream << std::endl << "Owned Regular Vertex List";
debugStream << std::endl << "Owned Regular Vertex List" << std::endl;
vtkm::worklet::contourtree_augmented::PrintHeader(globalIds.GetNumberOfValues(), debugStream);
vtkm::worklet::contourtree_augmented::PrintIndices("GlobalId", globalIds, debugStream);
vtkm::worklet::contourtree_augmented::PrintIndices("GlobalId", globalIds, -1, debugStream);
VTKM_LOG_S(vtkm::cont::LogLevel::Info, debugStream.str());
}
#endif
@ -268,19 +265,21 @@ void HierarchicalHyperSweeper<MeshType, FieldType>::InitializeIntrinsicVertexCou
auto findRegularByGlobal = hierarchicalTree.GetFindRegularByGlobal();
auto computeSuperparentIdsWorklet = vtkm::worklet::contourtree_distributed::
hierarchical_hyper_sweeper::InitializeIntrinsicVertexCountComputeSuperparentIdsWorklet();
localInvoke(computeSuperparentIdsWorklet, // worklet to run
globalIds, // input
findRegularByGlobal, // input
hierarchicalTree.Regular2Supernode, // input
hierarchicalTree.Superparents, // input
superparents // output
std::cout << "Invoking computeSuperparentIdsWorklet" << std::endl;
Invoke(computeSuperparentIdsWorklet, // worklet to run
globalIds, // input
findRegularByGlobal, // input
hierarchicalTree.Regular2Supernode, // input
hierarchicalTree.Superparents, // input
superparents // output
);
}
#ifdef DEBUG_PRINT
{
std::stringstream debugStream;
vtkm::worklet::contourtree_augmented::PrintIndices("Superparents", superparents, debugStream);
vtkm::worklet::contourtree_augmented::PrintIndices(
"Superparents", superparents, -1, debugStream);
VTKM_LOG_S(vtkm::cont::LogLevel::Info, debugStream.str());
}
#endif
@ -291,41 +290,34 @@ void HierarchicalHyperSweeper<MeshType, FieldType>::InitializeIntrinsicVertexCou
#ifdef DEBUG_PRINT
{
std::stringstream debugStream;
vtkm::worklet::contourtree_augmented::PrintIndices("Sorted SP", superparents, debugStream);
vtkm::worklet::contourtree_augmented::PrintIndices("Sorted SP", superparents, -1, debugStream);
VTKM_LOG_S(vtkm::cont::LogLevel::Info, debugStream.str());
}
#endif
// initialize the counts to zero.
vtkm::cont::Algorithm::Copy(
vtkm::cont::ArrayHandleConstant<vtkm::Id>(static_cast<vtkm::Id>(0),
hierarchicalTree.Supernodes.GetNumberOfValues()),
superarcRegularCounts);
{ // scope to make sure temporary variables are deleted
vtkm::worklet::contourtree_distributed::hierarchical_hyper_sweeper::
InitializeIntrinsicVertexCountInitalizeCountsWorklet initalizeCountsWorklet;
// set the count to the Id one off the high end of each range
localInvoke(initalizeCountsWorklet, // worklet
superparents, // input domain
superarcRegularCounts // output
);
}
vtkm::cont::Algorithm::Fill(
superarcRegularCounts, vtkm::Id{ 0 }, this->HierarchicalTree.Supernodes.GetNumberOfValues());
// set the count to the Id one off the high end of each range
Invoke(vtkm::worklet::contourtree_distributed::hierarchical_hyper_sweeper::
InitializeIntrinsicVertexCountInitalizeCountsWorklet{},
superparents, // input domain
superarcRegularCounts // output
);
// now repeat to subtract out the low end
{
vtkm::worklet::contourtree_distributed::hierarchical_hyper_sweeper::
InitializeIntrinsicVertexCountSubtractLowEndWorklet subtractLowEndWorklet;
localInvoke(subtractLowEndWorklet, // worklet
superparents, // input domain
superarcRegularCounts // output
);
}
Invoke(vtkm::worklet::contourtree_distributed::hierarchical_hyper_sweeper::
InitializeIntrinsicVertexCountSubtractLowEndWorklet{},
superparents, // input domain
superarcRegularCounts // output
);
// and that is that
#ifdef DEBUG_PRINT
{
std::stringstream debugStream;
vtkm::worklet::contourtree_augmented::PrintIndices(
"SuperarcRegularCounts", superarcRegularCounts, debugStream);
"SuperarcRegularCounts", superarcRegularCounts, -1, debugStream);
VTKM_LOG_S(vtkm::cont::LogLevel::Info, debugStream.str());
}
#endif
@ -333,36 +325,38 @@ void HierarchicalHyperSweeper<MeshType, FieldType>::InitializeIntrinsicVertexCou
// routine to do the local hypersweep using addition / subtraction
template <typename MeshType, typename FieldType>
void HierarchicalHyperSweeper<MeshType, FieldType>::LocalHyperSweep()
template <typename SweepValueType, typename ContourTreeFieldType>
void HierarchicalHyperSweeper<SweepValueType, ContourTreeFieldType>::LocalHyperSweep()
{ // LocalHyperSweep()
// TODO: Implement this function
#ifdef DEBUG_PRINT
VTKM_LOG_S(vtkm::cont::LogLevel::Info,
DebugPrint(std::string("Hypersweep Block ") + std::to_string(blockId) +
DebugPrint(std::string("Hypersweep Block ") + std::to_string(BlockId) +
std::string(" Starting Local HyperSweep"),
__FILE__,
__LINE__));
#endif
// I. Iterate over all rounds of the hyperstructure
for (vtkm::Id round = 0; round <= this->hierarchicalTree.NumRounds; round++)
for (vtkm::Id round = 0; round <= this->HierarchicalTree.NumRounds; round++)
{ // per round
#ifdef DEBUG_PRINT
VTKM_LOG_S(vtkm::cont::LogLevel::Info,
DebugPrint(std::string("Hypersweep Block ") + std::to_string(blockId) +
DebugPrint(std::string("Hypersweep Block ") + std::to_string(BlockId) +
std::string(" Round ") + std::to_string(round) +
std::string(" Step 0 Starting Round"),
__FILE__,
__LINE__));
#endif
// A. Iterate over all iterations of the round
for (vtkm::Id iteration = 0; iteration < this->HierarchicalTree.NumIterations[round];
iteration++)
auto numIterationsPortal =
this->HierarchicalTree.NumIterations
.ReadPortal(); // TODO/FIXME: Use portal? Or something more efficient?
for (vtkm::Id iteration = 0; iteration < numIterationsPortal.Get(round); iteration++)
{ // per iteration
#ifdef DEBUG_PRINT
VTKM_LOG_S(vtkm::cont::LogLevel::Info,
DebugPrint(std::string("Hypersweep Block ") + std::to_string(blockId) +
DebugPrint(std::string("Hypersweep Block ") + std::to_string(BlockId) +
std::string(" Round ") + std::to_string(round) +
std::string(" Step 1 Iteration ") + std::to_string(iteration) +
std::string(" Step A Starting Iteration"),
@ -370,16 +364,18 @@ void HierarchicalHyperSweeper<MeshType, FieldType>::LocalHyperSweep()
__LINE__));
#endif
// 1. Establish the range of supernode Ids that we want to process
vtkm::Id firstSupernode = this->HierarchicalTree.FirstSupernodePerIteration[round][iteration];
vtkm::Id lastSupernode =
this->HierarchicalTree.FirstSupernodePerIteration[round][iteration + 1];
// TODO/FIXME: Use portal? Or is there a more efficient way?
auto firstSupernodePerIterationPortal =
this->HierarchicalTree.FirstSupernodePerIteration[round].ReadPortal();
vtkm::Id firstSupernode = firstSupernodePerIterationPortal.Get(iteration);
vtkm::Id lastSupernode = firstSupernodePerIterationPortal.Get(iteration + 1);
// call the routine that computes the dependent weights for each superarc in that range
this->ComputeSuperarcDependentWeights(round, iteration, firstSupernode, lastSupernode);
#ifdef DEBUG_PRINT
VTKM_LOG_S(vtkm::cont::LogLevel::Info,
DebugPrint(std::string("Hypersweep Block ") + std::to_string(blockId) +
DebugPrint(std::string("Hypersweep Block ") + std::to_string(BlockId) +
std::string(" Round ") + std::to_string(round) +
std::string(" Step 1 Iteration ") + std::to_string(iteration) +
std::string(" Step B Dependent Weights Computed"),
@ -391,7 +387,7 @@ void HierarchicalHyperSweeper<MeshType, FieldType>::LocalHyperSweep()
#ifdef DEBUG_PRINT
VTKM_LOG_S(vtkm::cont::LogLevel::Info,
DebugPrint(std::string("Hypersweep Block ") + std::to_string(blockId) +
DebugPrint(std::string("Hypersweep Block ") + std::to_string(BlockId) +
std::string(" Round ") + std::to_string(round) +
std::string(" Step 1 Iteration ") + std::to_string(iteration) +
std::string(" Step C Transfer Weights Computed"),
@ -404,7 +400,7 @@ void HierarchicalHyperSweeper<MeshType, FieldType>::LocalHyperSweep()
#ifdef DEBUG_PRINT
VTKM_LOG_S(vtkm::cont::LogLevel::Info,
DebugPrint(std::string("Hypersweep Block ") + std::to_string(blockId) +
DebugPrint(std::string("Hypersweep Block ") + std::to_string(BlockId) +
std::string(" Round ") + std::to_string(round) +
std::string(" Step 1 Iteration ") + std::to_string(iteration) +
std::string(" Step D Weights Transferred"),
@ -415,7 +411,7 @@ void HierarchicalHyperSweeper<MeshType, FieldType>::LocalHyperSweep()
#ifdef DEBUG_PRINT
VTKM_LOG_S(vtkm::cont::LogLevel::Info,
DebugPrint(std::string("Hypersweep Block ") + std::to_string(blockId) +
DebugPrint(std::string("Hypersweep Block ") + std::to_string(BlockId) +
std::string(" Round ") + std::to_string(round) +
std::string(" Step 2 Ending Round"),
__FILE__,
@ -426,30 +422,29 @@ void HierarchicalHyperSweeper<MeshType, FieldType>::LocalHyperSweep()
// routine to compute the correct weights dependent on each superarc in a subrange (defined by the round & iteration)
template <typename MeshType, typename FieldType>
void HierarchicalHyperSweeper<MeshType, FieldType>::ComputeSuperarcDependentWeights(
vtkm::Id round,
vtkm::Id iteration,
vtkm::Id firstSupernode,
vtkm::Id lastSupernode)
template <typename SweepValueType, typename ContourTreeFieldType>
void HierarchicalHyperSweeper<SweepValueType, ContourTreeFieldType>::
ComputeSuperarcDependentWeights(
vtkm::Id round,
vtkm::Id, // iteration, // Kept parameter in case we need it for debugging.
vtkm::Id firstSupernode,
vtkm::Id lastSupernode)
{ // ComputeSuperarcDependentWeights()
(void)
iteration; // avoid compiler warning for unused parmeter. Kept parameter in case we need it for debugging.
vtkm::Id numSupernodesToProcess = lastSupernode - firstSupernode;
// 2. Use sorted prefix sum to compute the total weight to contribute to the super/hypertarget
// Same as std::partial_sum(sweepValues.begin() + firstSupernode, sweepValues.begin() + lastSupernode, valuePrefixSum.begin() + firstSupernode);
{
vtkm::Id numValuesToCopy = lastSupernode - firstSupernode;
// DependentValues[firstSuperNode, lastSupernode)
vtkm::cont::ArrayHandleView<vtkm::worklet::contourtree_augmented::IdArrayType>
dependentValuesView(this->DependentValues, // subset DependentValues
firstSupernode, // start at firstSupernode
numValuesToCopy); // until lastSuperNode (not inclued)
dependentValuesView(this->DependentValues, // subset DependentValues
firstSupernode, // start at firstSupernode
numSupernodesToProcess); // until lastSuperNode (not inclued)
// Target array
vtkm::cont::ArrayHandleView<vtkm::worklet::contourtree_augmented::IdArrayType>
valuePrefixSumView(this->ValuePrefixSum, // subset ValuePrefixSum
firstSupernode, // start at firstSupernode
numValuesToCopy); // until lastSuperNode (not inclued)
valuePrefixSumView(this->ValuePrefixSum, // subset ValuePrefixSum
firstSupernode, // start at firstSupernode
numSupernodesToProcess); // until lastSuperNode (not inclued)
// Compute the partial sum for DependentValues[firstSuperNode, lastSupernode) and write to ValuePrefixSum[firstSuperNode, lastSupernode)
vtkm::cont::Algorithm::ScanInclusive(dependentValuesView, // input
valuePrefixSumView); // result of partial sum
@ -466,23 +461,22 @@ void HierarchicalHyperSweeper<MeshType, FieldType>::ComputeSuperarcDependentWeig
// 3. Compute the segmented weights from the prefix sum array
{
// Create views of the subranges of the arrays we need to update
vtkm::Id numValuesToProcess = lastSupernode - firstSupernode;
vtkm::cont::ArrayHandleCounting<vtkm::Id> supernodeIndex(
firstSupernode, static_cast<vtkm::Id>(1), numValuesToProcess);
firstSupernode, vtkm::Id{ 1 }, numSupernodesToProcess);
vtkm::cont::ArrayHandleView<vtkm::worklet::contourtree_augmented::IdArrayType>
hierarchicalTreeSuperarcsView(
this->HierarchicalTree.Superarcs, firstSupernode, numValuesToProcess);
this->HierarchicalTree.Superarcs, firstSupernode, numSupernodesToProcess);
vtkm::cont::ArrayHandleView<vtkm::worklet::contourtree_augmented::IdArrayType>
hierarchicalTreeHyperparentsView(
this->HierarchicalTree.Hyperparents, firstSupernode, numValuesToProcess);
this->HierarchicalTree.Hyperparents, firstSupernode, numSupernodesToProcess);
vtkm::cont::ArrayHandleView<vtkm::worklet::contourtree_augmented::IdArrayType>
hierarchicalTreeHypernodesView(
this->HierarchicalTree.Hypernodes, firstSupernode, numValuesToProcess);
this->HierarchicalTree.Hypernodes, firstSupernode, numSupernodesToProcess);
vtkm::cont::ArrayHandleView<vtkm::worklet::contourtree_augmented::IdArrayType>
dependentValuesView(this->DependentValues, firstSupernode, numValuesToProcess);
dependentValuesView(this->DependentValues, firstSupernode, numSupernodesToProcess);
// create the worklet
vtkm::worklet::contourtree_distributed::hierarchical_hyper_sweeper::
ComputeSuperarcDependentWeightsWorklet<FieldType>
ComputeSuperarcDependentWeightsWorklet<SweepValueType>
computeSuperarcDependentWeightsWorklet(
firstSupernode, round, this->HierarchicalTree.NumRounds);
// Execute the worklet
@ -491,7 +485,7 @@ void HierarchicalHyperSweeper<MeshType, FieldType>::ComputeSuperarcDependentWeig
supernodeIndex, // input counting index [firstSupernode, lastSupernode)
hierarchicalTreeSuperarcsView, // input view of hierarchicalTree.Superarcs[firstSupernode, lastSupernode)
hierarchicalTreeHyperparentsView, // input view of hierarchicalTree.Hyperparents[firstSupernode, lastSupernode)
this->hierarchicalTree.Hypernodes, // input full hierarchicalTree.Hypernodes array
this->HierarchicalTree.Hypernodes, // input full hierarchicalTree.Hypernodes array
this->ValuePrefixSum, // input full ValuePrefixSum array
dependentValuesView // output view of sweepValues[firstSu
);
@ -500,35 +494,32 @@ void HierarchicalHyperSweeper<MeshType, FieldType>::ComputeSuperarcDependentWeig
// routine to compute the weights to transfer to superarcs (defined by the round & iteration)
template <typename MeshType, typename FieldType>
void HierarchicalHyperSweeper<MeshType, FieldType>::ComputeSuperarcTransferWeights(
template <typename SweepValueType, typename ContourTreeFieldType>
void HierarchicalHyperSweeper<SweepValueType, ContourTreeFieldType>::ComputeSuperarcTransferWeights(
vtkm::Id round,
vtkm::Id iteration,
vtkm::Id, // iteration, // Kept parameter in case we need it for debugging.
vtkm::Id firstSupernode,
vtkm::Id lastSupernode)
{ // ComputeSuperarcTransferWeights()
(void)
iteration; // avoid compiler warning for unused parmeter. Kept parameter in case we need it for debugging.
// At this stage, we would otherwise transfer weights by hyperarc, but attachment points don't *have* hyperarcs
// so we will do a transfer by superarc instead, making sure that we only transfer from the last superarc in each
// hyperarc, plus for any attachment point
vtkm::Id numSupernodesToProcess = lastSupernode - firstSupernode;
// 4. Set the amount each superarc wants to transfer, reusing the valuePrefixSum array for the purpose
// and the transfer target
{ // scope ComputeSuperarcTransferWeightsWorklet to make sure temp variables are cleared
// Create ArrayHandleViews of the subrange of values that we need to update
vtkm::Id numValuesToProcess = lastSupernode - firstSupernode;
vtkm::cont::ArrayHandleCounting<vtkm::Id> supernodeIndex(
firstSupernode, static_cast<vtkm::Id>(1), numValuesToProcess);
firstSupernode, vtkm::Id{ 1 }, numSupernodesToProcess);
vtkm::cont::ArrayHandleView<vtkm::worklet::contourtree_augmented::IdArrayType>
hierarchicalTreeSupernodesView(
this->HierarchicalTree.Supernodes, firstSupernode, numValuesToProcess);
this->HierarchicalTree.Supernodes, firstSupernode, numSupernodesToProcess);
vtkm::cont::ArrayHandleView<vtkm::worklet::contourtree_augmented::IdArrayType>
hierarchicalTreeSuperarcsView(
this->HierarchicalTree.Superarcs, firstSupernode, numValuesToProcess);
this->HierarchicalTree.Superarcs, firstSupernode, numSupernodesToProcess);
vtkm::cont::ArrayHandleView<vtkm::worklet::contourtree_augmented::IdArrayType>
transferTargetView(this->transferTarget, firstSupernode, numValuesToProcess);
transferTargetView(this->TransferTarget, firstSupernode, numSupernodesToProcess);
// instantiate the worklet
vtkm::worklet::contourtree_distributed::hierarchical_hyper_sweeper::
ComputeSuperarcTransferWeightsWorklet computeSuperarcTransferWeightsWorklet(
@ -548,12 +539,11 @@ void HierarchicalHyperSweeper<MeshType, FieldType>::ComputeSuperarcTransferWeigh
// 5. Now we need to sort the transfer targets into contiguous segments
{
// create view of superSortPermute[firstSupernode, lastSupernode) for sorting
vtkm::Id numValuesToProcess = lastSupernode - firstSupernode;
vtkm::cont::ArrayHandleView<vtkm::worklet::contourtree_augmented::IdArrayType>
superSortPermuteView(this->SuperSortPermute, firstSupernode, numValuesToProcess);
superSortPermuteView(this->SuperSortPermute, firstSupernode, numSupernodesToProcess);
// create comperator for the sort
vtkm::worklet::contourtree_distributed::hierarchical_hyper_sweeper::TransferTargetComperator
transferTargetComperator(this->transferTarget);
transferTargetComperator(this->TransferTarget);
// sort the subrange of our array
vtkm::cont::Algorithm::Sort(superSortPermuteView, transferTargetComperator);
}
@ -567,11 +557,10 @@ void HierarchicalHyperSweeper<MeshType, FieldType>::ComputeSuperarcTransferWeigh
// }
{
// copy transfer target in the sorted order
vtkm::Id numValuesToProcess = lastSupernode - firstSupernode;
vtkm::cont::ArrayHandleView<vtkm::worklet::contourtree_augmented::IdArrayType>
sortedTransferTargetView(this->SortedTransferTarget, firstSupernode, numValuesToProcess);
sortedTransferTargetView(this->SortedTransferTarget, firstSupernode, numSupernodesToProcess);
vtkm::cont::ArrayHandleView<vtkm::worklet::contourtree_augmented::IdArrayType>
superSortPermuteView(this->SuperSortPermute, firstSupernode, numValuesToProcess);
superSortPermuteView(this->SuperSortPermute, firstSupernode, numSupernodesToProcess);
auto permutedTransferTarget =
vtkm::cont::make_ArrayHandlePermutation(superSortPermuteView, // idArray
this->TransferTarget); // valueArray
@ -579,7 +568,7 @@ void HierarchicalHyperSweeper<MeshType, FieldType>::ComputeSuperarcTransferWeigh
// Note that any values associated with NO_SUCH_ELEMENT will be ignored
// copy transfer weight in the sorted order
vtkm::cont::ArrayHandleView<vtkm::worklet::contourtree_augmented::IdArrayType>
valuePrefixSumView(this->ValuePrefixSum, firstSupernode, numValuesToProcess);
valuePrefixSumView(this->ValuePrefixSum, firstSupernode, numSupernodesToProcess);
auto permutedDependentValues =
vtkm::cont::make_ArrayHandlePermutation(superSortPermuteView, // idArray
this->DependentValues); // valueArray
@ -589,31 +578,28 @@ void HierarchicalHyperSweeper<MeshType, FieldType>::ComputeSuperarcTransferWeigh
// routine to transfer the weights
template <typename MeshType, typename FieldType>
void HierarchicalHyperSweeper<MeshType, FieldType>::TransferWeights(vtkm::Id round,
vtkm::Id iteration,
vtkm::Id firstSupernode,
vtkm::Id lastSupernode)
template <typename SweepValueType, typename ContourTreeFieldType>
void HierarchicalHyperSweeper<SweepValueType, ContourTreeFieldType>::TransferWeights(
vtkm::Id, // round, // Kept parameters in case we need it for debugging.
vtkm::Id, // iteration, // Kept parameters in case we need it for debugging.
vtkm::Id firstSupernode,
vtkm::Id lastSupernode)
{ // TransferWeights()
// avoid compiler warning for unused parmeters. Kept parameters in case we need it for debugging.
(void)round;
(void)iteration;
// 7. Now perform a segmented prefix sum
vtkm::Id numSupernodesToProcess = lastSupernode - firstSupernode;
// Same as std::partial_sum(valuePrefixSum.begin() + firstSupernode, valuePrefixSum.begin() + lastSupernode, valuePrefixSum.begin() + firstSupernode);
{
vtkm::Id numValuesToCopy = lastSupernode - firstSupernode;
// ValuePrefixSum[firstSuperNode, lastSupernode)
vtkm::cont::ArrayHandleView<vtkm::worklet::contourtree_augmented::IdArrayType>
valuePrefixSumView(this->ValuePrefixSum, // subset ValuePrefixSum
firstSupernode, // start at firstSupernode
numValuesToCopy); // until lastSuperNode (not inclued)
valuePrefixSumView(this->ValuePrefixSum, // subset ValuePrefixSum
firstSupernode, // start at firstSupernode
numSupernodesToProcess); // until lastSuperNode (not inclued)
// TODO: If it is safe to use the same array as input and output for ScanInclusive then this code should be updated to avoid the extra copy
// In this case our traget array is the same as our source array. For safety we
// store the values of our prefix sum in a temporary arrya and then copy the values
// back into our valuePrefixSumView at the end
vtkm::cont::ArrayHandle<FieldType> tempScanInclusiveTarget;
tempScanInclusiveTarget.Allocate(numValuesToCopy);
vtkm::worklet::contourtree_augmented::IdArrayType tempScanInclusiveTarget;
tempScanInclusiveTarget.Allocate(numSupernodesToProcess);
// Compute the partial sum for DependentValues[firstSuperNode, lastSupernode) and write to ValuePrefixSum[firstSuperNode, lastSupernode)
vtkm::cont::Algorithm::ScanInclusive(valuePrefixSumView, // input
tempScanInclusiveTarget); // result of partial sum
@ -623,60 +609,54 @@ void HierarchicalHyperSweeper<MeshType, FieldType>::TransferWeights(vtkm::Id rou
// 7a. and 7b.
{
// Prepare the approbriate array views for our worklet. This is done to allow us to
// use FieldIn instead of having to transfer the entire array to the device when we
// really only need a subrange
vtkm::Id numValuesToProcess = lastSupernode - firstSupernode;
vtkm::cont::ArrayHandleCounting<vtkm::Id> supernodeIndex(
firstSupernode, static_cast<vtkm::Id>(1), numValuesToProcess);
vtkm::cont::ArrayHandleView<vtkm::worklet::contourtree_augmented::IdArrayType>
sortedTransferTargetView(this->SortedTransferTarget, firstSupernode, numValuesToProcess);
vtkm::cont::ArrayHandleView<vtkm::worklet::contourtree_augmented::IdArrayType>
sortedTransferTargetShiftedView(
this->SortedTransferTarget, firstSupernode + 1, numValuesToProcess);
vtkm::cont::ArrayHandleView<vtkm::worklet::contourtree_augmented::IdArrayType>
valuePrefixSumView(this->ValuePrefixSum, firstSupernode, numValuesToProcess);
vtkm::cont::ArrayHandleView<vtkm::worklet::contourtree_augmented::IdArrayType>
valuePrefixSumShiftedView(this->ValuePrefixSum, firstSupernode - 1, numValuesToProcess);
auto sweepValuePermuted =
vtkm::cont::make_ArrayHandlePermutation(sortedTransferTargetView, // idArray
this->DependentValues); // valueArray
// 7a. Find the RHE of each group and transfer the prefix sum weight
// Note that we do not compute the transfer weight separately, we add it in place instead
// Instantiate the worklet
auto supernodeIndex =
vtkm::cont::make_ArrayHandleCounting(firstSupernode, vtkm::Id{ 1 }, numSupernodesToProcess);
VTKM_ASSERT(firstSupernode + numSupernodesToProcess <=
this->ValuePrefixSum.GetNumberOfValues());
auto valuePrefixSumView = vtkm::cont::make_ArrayHandleView(
this->ValuePrefixSum, firstSupernode, numSupernodesToProcess);
vtkm::worklet::contourtree_distributed::hierarchical_hyper_sweeper::
TransferWeightsUpdateRHEWorklet transferWeightsUpdateRHEWorklet(lastSupernode);
// Invoke the worklet
this->Invoke(
transferWeightsUpdateRHEWorklet, // worklet
supernodeIndex, // input counting array [firstSupernode, lastSupernode)
sortedTransferTargetView, // input view of sortedTransferTarget[firstSupernode, lastSupernode)
sortedTransferTargetShiftedView, // input view of sortedTransferTarget[firstSupernode+1, lastSupernode+1)
valuePrefixSumView, // input view of valuePrefixSum[firstSupernode, lastSupernode)
sweepValuePermuted // output view of sweepValues permuted by sortedTransferTarget[firstSupernode, lastSupernode). Use FieldInOut since we don't overwrite all values.
);
this->Invoke(transferWeightsUpdateRHEWorklet, // worklet
supernodeIndex, // input counting array [firstSupernode, lastSupernode)
this->SortedTransferTarget,
valuePrefixSumView, // input view of valuePrefixSum[firstSupernode, lastSupernode)
this->DependentValues);
}
{
VTKM_ASSERT(firstSupernode + 1 + numSupernodesToProcess - 1 <=
this->SortedTransferTarget.GetNumberOfValues());
auto sortedTransferTargetView = vtkm::cont::make_ArrayHandleView(
this->SortedTransferTarget, firstSupernode + 1, numSupernodesToProcess - 1);
VTKM_ASSERT(firstSupernode + 1 + numSupernodesToProcess - 1 <=
this->SortedTransferTarget.GetNumberOfValues());
auto sortedTransferTargetShiftedView = vtkm::cont::make_ArrayHandleView(
this->SortedTransferTarget, firstSupernode + 1, numSupernodesToProcess - 1);
auto valuePrefixSumPreviousValueView = vtkm::cont::make_ArrayHandleView(
this->ValuePrefixSum, firstSupernode, numSupernodesToProcess - 1);
// 7b. Now find the LHE of each group and subtract out the prior weight
vtkm::worklet::contourtree_distributed::hierarchical_hyper_sweeper::
TransferWeightsUpdateLHEWorklet transferWeightsUpdateLHEWorklet(firstSupernode);
this->Invoke(
transferWeightsUpdateLHEWorklet, // worklet
supernodeIndex, // input counting array [firstSupernode, lastSupernode)
sortedTransferTargetView, // input view of sortedTransferTarget[firstSupernode, lastSupernode)
sortedTransferTargetShiftedView, // input view of sortedTransferTarget[firstSupernode+1, lastSupernode+1)
valuePrefixSumShiftedView, // input view of valuePrefixSum[firstSupernode-1, lastSupernode-1)
sweepValuePermuted // output view of sweepValues permuted by sortedTransferTarget[firstSupernode, lastSupernode). Use FieldInOut since we don't overwrite all values.
);
this->Invoke(vtkm::worklet::contourtree_distributed::hierarchical_hyper_sweeper::
TransferWeightsUpdateLHEWorklet{},
sortedTransferTargetView,
sortedTransferTargetShiftedView,
valuePrefixSumPreviousValueView,
this->DependentValues);
}
} // TransferWeights()
// debug routine
template <typename MeshType, typename FieldType>
std::string HierarchicalHyperSweeper<MeshType, FieldType>::DebugPrint(std::string message,
const char* fileName,
long lineNum) const
template <typename SweepValueType, typename ContourTreeFieldType>
std::string HierarchicalHyperSweeper<SweepValueType, ContourTreeFieldType>::DebugPrint(
std::string message,
const char* fileName,
long lineNum) const
{ // DebugPrint()
std::stringstream resultStream;
resultStream << std::endl;
@ -695,20 +675,20 @@ std::string HierarchicalHyperSweeper<MeshType, FieldType>::DebugPrint(std::strin
vtkm::worklet::contourtree_augmented::PrintIndices(
"Dependent", this->DependentValues, -1, resultStream);
vtkm::worklet::contourtree_augmented::PrintIndices(
"Prefix Sum", this->ValuePrefixSum - 1, resultStream);
"Prefix Sum", this->ValuePrefixSum, -1, resultStream);
vtkm::worklet::contourtree_augmented::PrintIndices(
"Transfer To", this->TransferTarget - 1, resultStream);
"Transfer To", this->TransferTarget, -1, resultStream);
vtkm::worklet::contourtree_augmented::PrintIndices(
"Sorted Transfer", this->SortedTransferTarget - 1, resultStream);
"Sorted Transfer", this->SortedTransferTarget, -1, resultStream);
vtkm::worklet::contourtree_augmented::PrintIndices(
"Sort Permute", this->SuperSortPermute - 1, resultStream);
"Sort Permute", this->SuperSortPermute, -1, resultStream);
return resultStream.str();
} // DebugPrint()
// Routine to save the hierarchical tree to file
template <typename MeshType, typename FieldType>
void HierarchicalHyperSweeper<MeshType, FieldType>::SaveHierarchicalContourTreeDot(
template <typename SweepValueType, typename ContourTreeFieldType>
void HierarchicalHyperSweeper<SweepValueType, ContourTreeFieldType>::SaveHierarchicalContourTreeDot(
std::string message,
const char* outFileName) const
{ // SaveHierarchicalContourTreeDot()

4
vtkm/worklet/contourtree_distributed/hierarchical_hyper_sweeper/ComputeSuperarcDependentWeightsWorklet.h
View File

@ -75,8 +75,6 @@ public:
WholeArrayIn valuePrefixSum, // whole valuePrefixSum array
FieldInOut dependentValuesView // output view of dependentValues[firstSupernode, lastSupernode)
);
using ExecutionSignature = void(InputIndex, _1, _2);
using InputDomain = _1;
// Default Constructor
VTKM_EXEC_CONT
@ -89,7 +87,7 @@ public:
{
}
template <typename InFieldPortalType, typename InOutFieldPortalType>
template <typename InFieldPortalType>
VTKM_EXEC void operator()(
const vtkm::Id& supernode,
const vtkm::Id& superarcTo, // same as hierarchicalTree.superarcs[supernode];

16
vtkm/worklet/contourtree_distributed/hierarchical_hyper_sweeper/ComputeSuperarcTransferWeightsWorklet.h
View File

@ -76,12 +76,10 @@ public:
hierarchicalTreeSupernodesView, // view of hierarchicalTree.supernodes[firstSupernode, lastSupernode)
WholeArrayIn hierarchicalTreeSuperparents, // whole array of hierarchicalTree.superparents
WholeArrayIn hierarchicalTreeHyperparents, // whole array of hierarchicalTree.hyperparents
FieldInOut
FieldIn
hierarchicalTreeSuperarcsView, // view of hierarchicalTree.superarcs[firstSupernode, lastSupernode)
FieldOut transferTargetView // view of transferTarget[firstSupernode, lastSupernode)
);
using ExecutionSignature = void(_1, _2, _3, _4, _5, _6);
using InputDomain = _1;
// Default Constructor
VTKM_EXEC_CONT
@ -94,14 +92,14 @@ public:
{
}
template <typename InFieldPortalType, typename InOutFieldPortalType>
template <typename InFieldPortalType>
VTKM_EXEC void operator()(
const vtkm::Id& supernode,
const vtkm::Id& supernodeRegularId, // same as hierarchicalTree.supernodes[supernode];
const InFieldPortalType& hierarchicalTreeSuperparentsPortal,
const InFieldPortalType& hierarchicalTreeHyperparentsPortal,
vtkm::Id& superarcTo, // same as hierarchicalTree.superarcs[supernode];
vtkm::Id& transferTarget // same as transferTarget[supernode]
const vtkm::Id& superarcTo, // same as hierarchicalTree.superarcs[supernode];
vtkm::Id& transferTarget // same as transferTarget[supernode]
) const
{
// per supernode
@ -131,10 +129,8 @@ public:
} // not a superarc we care about
else
{ // a superarc we care about
// strip off the flag bits
superarcTo = vtkm::worklet::contourtree_augmented::MaskedIndex(superarcTo);
// and set the weight & target
transferTarget = superarcTo;
// strip off the flag bits and set the weight & target
transferTarget = vtkm::worklet::contourtree_augmented::MaskedIndex(superarcTo);
} // a superarc we care about
} // actual superarc

10
vtkm/worklet/contourtree_distributed/hierarchical_hyper_sweeper/InitializeIntrinsicVertexCountComputeSuperparentIdsWorklet.h
View File

@ -87,11 +87,11 @@ public:
InitializeIntrinsicVertexCountComputeSuperparentIdsWorklet() {}
template <typename ExecObjType, typename InFieldPortalType>
VTKM_EXEC vtkm::Id operator()(const vtkm::Id& globalId,
const ExecObjType& findRegularByGlobal,
const InFieldPortalType& hierarchicalTreeRegular2SupernodePortal,
const InFieldPortalType& hierarchicalTreeSuperparentsPortal,
vtkm::Id& superparent) const
VTKM_EXEC void operator()(const vtkm::Id& globalId,
const ExecObjType& findRegularByGlobal,
const InFieldPortalType& hierarchicalTreeRegular2SupernodePortal,
const InFieldPortalType& hierarchicalTreeSuperparentsPortal,
vtkm::Id& superparent) const
{
// per vertex
// retrieve the regular Id (should ALWAYS exist)

23
vtkm/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferTargetComperator.h
View File

@ -55,6 +55,7 @@
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ExecutionObjectBase.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
namespace vtkm
{
@ -84,7 +85,27 @@ public:
VTKM_EXEC
bool operator()(const vtkm::Id& left, const vtkm::Id& right) const
{ // operator()
return this->SuperarcPortal.Get(left) < this->SuperarcPortal.Get(right);
// NOTE: We need to explicitly check for NO_SUCH_ELEMENT here since vtkm::Id is signed
// while the index time in PPP is unsigned. Thus, for PPP "regular" indices are always
// smaller that NO_SUCH_ELEMENT, while with the signed vtkm::Id, NO_SUCH_ELEMENT is
// negative and the order is not as intented.
// TODO/FIXME: Verify this implementation is correct.
// TODO/FIXME: Is there a better way to do this?
auto leftVal = this->SuperarcPortal.Get(left);
auto rightVal = this->SuperarcPortal.Get(right);
if (vtkm::worklet::contourtree_augmented::NoSuchElement(leftVal))
{
return false;
}
else if (vtkm::worklet::contourtree_augmented::NoSuchElement(rightVal))
{
return true;
}
else
{
return vtkm::worklet::contourtree_augmented::MaskedIndex(leftVal) <
vtkm::worklet::contourtree_augmented::MaskedIndex(rightVal);
}
} // operator()
private:

70
vtkm/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateLHEWorklet.h
View File

@ -70,56 +70,32 @@ namespace hierarchical_hyper_sweeper
class TransferWeightsUpdateLHEWorklet : public vtkm::worklet::WorkletMapField
{
public:
using ControlSignature = void(
FieldIn supernodeIndex, // input counting array [firstSupernode, lastSupernode)
FieldIn
sortedTransferTargetView, // input view of sortedTransferTarget[firstSupernode, lastSupernode)
FieldIn
sortedTransferTargetShiftedView, // input view of sortedTransferTarget[firstSupernode+1, lastSupernode+1)
FieldIn
valuePrefixSumShiftedView, // input view of valuePrefixSum[firstSupernode-1, lastSupernode-1)
FieldInOut
dependentValuePermuted // output view of dependentValues permuted by sortedTransferTarget[firstSupernode, lastSupernode). Use FieldInOut since we don't overwrite all values.
);
using ExecutionSignature = void(_1, _2, _3, _4, _5);
using InputDomain = _1;
using ControlSignature = void(FieldIn sortedTransferTargetPortal,
FieldIn sortedTransferTargetShiftedView,
FieldIn valuePrefixSumShiftedView,
WholeArrayInOut dependentValuesPortal);
// Default Constructor
VTKM_EXEC_CONT
TransferWeightsUpdateLHEWorklet(const vtkm::Id& firstSupernode)
: FirstSupernode(firstSupernode)
{
}
VTKM_EXEC void operator()(
const vtkm::Id& supernode,
const vtkm::Id& sortedTransferTargetValue, // same as sortedTransferTarget[supernode]
const vtkm::Id& sortedTransferTargetNextValue, // same as sortedTransferTarget[supernode+1]
const vtkm::Id& valuePrefixSumPreviousValue, // same as valuePrefixSum[supernode-1]
vtkm::Id& dependentValue // same as dependentValues[sortedTransferTarget[supernode]]
) const
template <typename InOutPortalType>
VTKM_EXEC void operator()(const vtkm::Id& sortedTransferTargetValue,
const vtkm::Id& sortedTransferTargetPreviusValue,
const vtkm::Id& valuePrefixSumPreviousValue,
InOutPortalType& dependentValuesPortal) const
{
// per supernode
// ignore any that point at NO_SUCH_ELEMENT
if (vtkm::worklet::contourtree_augmented::NoSuchElement(sortedTransferTargetValue))
if (!vtkm::worklet::contourtree_augmented::NoSuchElement(sortedTransferTargetValue))
{
return;
if (sortedTransferTargetValue != sortedTransferTargetPreviusValue)
{
auto originalValue = dependentValuesPortal.Get(sortedTransferTargetValue);
dependentValuesPortal.Set(sortedTransferTargetValue,
originalValue - valuePrefixSumPreviousValue);
}
}
// the LHE at 0 is special - it subtracts zero. In practice, since NO_SUCH_ELEMENT will sort low, this will never
// occur, but let's keep the logic strict
if (supernode == this->FirstSupernode)
{ // LHE 0
dependentValue -= 0;
} // LHE 0
else if (sortedTransferTargetValue != sortedTransferTargetNextValue)
{ // LHE not 0
dependentValue -= valuePrefixSumPreviousValue;
} // LHE not 0
// In serial this worklet implements the following operation
/*
for (vtkm::Id supernode = firstSupernode; supernode < lastSupernode; supernode++)
for (vtkm::Id supernode = firstSupernode + 1; supernode < lastSupernode; supernode++)
{ // per supernode
// ignore any that point at NO_SUCH_ELEMENT
if (noSuchElement(sortedTransferTarget[supernode]))
@ -127,22 +103,14 @@ public:
// the LHE at 0 is special - it subtracts zero. In practice, since NO_SUCH_ELEMENT will sort low, this will never
// occur, but let's keep the logic strict
if (supernode == firstSupernode)
{ // LHE 0
dependentValues[sortedTransferTarget[supernode]] -= 0;
} // LHE 0
else if (sortedTransferTarget[supernode] != sortedTransferTarget[supernode-1])
if (sortedTransferTarget[supernode] != sortedTransferTarget[supernode-1])
{ // LHE not 0
dependentValues[sortedTransferTarget[supernode]] -= valuePrefixSum[supernode-1];
} // LHE not 0
} // per supernode
*/
} // operator()()
private:
const vtkm::Id FirstSupernode;
}; // TransferWeightsUpdateLHEWorklet
}; // TransferWeightsUpdateLHEWorklet
} // namespace hierarchical_hyper_sweeper
} // namespace contourtree_distributed

46
vtkm/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateRHEWorklet.h
View File

@ -71,18 +71,11 @@ namespace hierarchical_hyper_sweeper
class TransferWeightsUpdateRHEWorklet : public vtkm::worklet::WorkletMapField
{
public:
using ControlSignature = void(
FieldIn supernodeIndex, // input counting array [firstSupernode, lastSupernode)
FieldIn
sortedTransferTargetView, // input view of sortedTransferTarget[firstSupernode, lastSupernode)
FieldIn
sortedTransferTargetShiftedView, // input view of sortedTransferTarget[firstSupernode+1, lastSupernode+1)
FieldIn valuePrefixSumView, // input view of valuePrefixSum[firstSupernode, lastSupernode)
FieldInOut
dependentValuePermuted // output view of dependentValues permuted by sortedTransferTarget[firstSupernode, lastSupernode). Use FieldInOut since we don't overwrite all values.
);
using ExecutionSignature = void(_1, _2, _3, _4, _5);
using InputDomain = _1;
using ControlSignature =
void(FieldIn supernodeIndex, // input counting array [firstSupernode, lastSupernode)
WholeArrayIn sortedTransferTarget,
FieldIn valuePrefixSumView, // input view of valuePrefixSum[firstSupernode, lastSupernode)
WholeArrayInOut dependentValuesPortal);
// Default Constructor
VTKM_EXEC_CONT
@ -91,26 +84,25 @@ public:
{
}
VTKM_EXEC void operator()(
const vtkm::Id& supernode,
const vtkm::Id& sortedTransferTargetValue, // same as sortedTransferTarget[supernode]
const vtkm::Id& sortedTransferTargetNextValue, // same as sortedTransferTarget[supernode+1]
const vtkm::Id& valuePrefixSum, // same as valuePrefixSum[supernode]
vtkm::Id& dependentValue // same as dependentValues[sortedTransferTarget[supernode]]
) const
template <typename InPortalType, typename OutPortalType>
VTKM_EXEC void operator()(const vtkm::Id& supernode,
const InPortalType& sortedTransferTargetPortal,
const vtkm::Id& valuePrefixSum, // same as valuePrefixSum[supernode]
OutPortalType& dependentValuesPortal) const
{
// per supernode
// ignore any that point at NO_SUCH_ELEMENT
if (vtkm::worklet::contourtree_augmented::NoSuchElement(sortedTransferTargetValue))
vtkm::Id transferTarget = sortedTransferTargetPortal.Get(supernode);
if (!vtkm::worklet::contourtree_augmented::NoSuchElement(transferTarget))
{
return;
// the RHE of each segment transfers its weight (including all irrelevant prefixes)
if ((supernode == this->LastSupernode - 1) ||
(transferTarget != sortedTransferTargetPortal.Get(supernode + 1)))
{ // RHE of segment
auto originalValue = dependentValuesPortal.Get(transferTarget);
dependentValuesPortal.Set(transferTarget, originalValue + valuePrefixSum);
} // RHE of segment
}
// the RHE of each segment transfers its weight (including all irrelevant prefixes)
if ((supernode == this->LastSupernode - 1) ||
(sortedTransferTargetValue != sortedTransferTargetNextValue))
{ // RHE of segment
dependentValue += valuePrefixSum;
} // RHE of segment
// In serial this worklet implements the following operation
/*

336
vtkm/worklet/testing/UnitTestContourTreeUniformDistributed.cxx
View File

@ -50,10 +50,23 @@
// Oliver Ruebel (LBNL)
//==============================================================================
#define DEBUG_PRINT
#include <vtkm/cont/testing/MakeTestDataSet.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/worklet/contourtree_augmented/DataSetMesh.h>
#include <vtkm/worklet/contourtree_augmented/PrintVectors.h>
#include <vtkm/worklet/contourtree_augmented/meshtypes/ContourTreeMesh.h>
#include <vtkm/worklet/contourtree_distributed/HierarchicalContourTree.h>
#include <vtkm/worklet/contourtree_distributed/HierarchicalHyperSweeper.h>
#include <vtkm/worklet/contourtree_distributed/SpatialDecomposition.h>
// clang-format off
VTKM_THIRDPARTY_PRE_INCLUDE
#include <vtkm/thirdparty/diy/diy.h>
VTKM_THIRDPARTY_POST_INCLUDE
// clang-format on
namespace
{
@ -85,13 +98,336 @@ void TestContourTreeMeshCombine(const std::string& mesh1_filename,
VTKM_TEST_ASSERT(contourTreeMesh2.MaxNeighbors == combinedContourTreeMesh.MaxNeighbors);
}
/*
template <typename PortalType, typename T>
static inline VTKM_CONT bool test_equal_portal_stl_vector(const PortalType1& portal,
const T[] array,
{
if (portal.GetNumberOfValues() != size)
{
return false;
}
for (vtkm::Id index = 0; index < portal.GetNumberOfValues(); index++)
{
if (!test_equal(portal.Get(index), array[index]))
{
return false;
}
}
return true;
}
*/
template <typename ContourTreeDataFieldType>
struct HyperSweepBlock
{
HyperSweepBlock(
const vtkm::Id& blockNo,
const vtkm::Id3& origin,
const vtkm::Id3& size,
const vtkm::Id3& globalSize,
const vtkm::worklet::contourtree_distributed::HierarchicalContourTree<ContourTreeDataFieldType>&
hierarchicalContourTree)
: BlockNo(blockNo)
, Origin(origin)
, Size(size)
, GlobalSize(globalSize)
, HierarchicalContourTree(hierarchicalContourTree)
{
}
// Mesh information
vtkm::Id BlockNo;
vtkm::Id3 Origin;
vtkm::Id3 Size;
vtkm::Id3 GlobalSize;
// Hierarchical contour tree for this block
const vtkm::worklet::contourtree_distributed::HierarchicalContourTree<ContourTreeDataFieldType>&
HierarchicalContourTree;
// Computed values
vtkm::cont::ArrayHandle<vtkm::Id> IntrinsicVolume;
vtkm::cont::ArrayHandle<vtkm::Id> DependentVolume;
};
template <typename ContourTreeDataFieldType>
struct CobmineHyperSweepBlockFunctor
{
void operator()(HyperSweepBlock<ContourTreeDataFieldType>* b,
const vtkmdiy::ReduceProxy& rp, // communication proxy
const vtkmdiy::RegularSwapPartners& // partners of the current block (unused)
) const
{
// Get our rank and DIY id
//const vtkm::Id rank = vtkm::cont::EnvironmentTracker::GetCommunicator().rank();
const auto selfid = rp.gid();
std::vector<int> incoming;
rp.incoming(incoming);
for (const int ingid : incoming)
{
std::cout << rp.round() << std::endl;
auto roundNo = rp.round() - 1;
// NOTE/IMPORTANT: In each round we should have only one swap partner (despite for-loop here).
// If that assumption does not hold, it will break things.
// NOTE/IMPORTANT: This assumption only holds if the number of blocks is a power of two.
// Otherwise, we may need to process more than one incoming block
if (ingid != selfid)
{
vtkm::Id incomingBlockNo;
rp.dequeue(ingid, incomingBlockNo);
// std::cout << "Combining block " << b->BlockNo << " with " << incomingBlockNo << std::endl;
vtkm::cont::ArrayHandle<vtkm::Id> incomingIntrinsicVolume;
rp.dequeue(ingid, incomingIntrinsicVolume);
vtkm::cont::ArrayHandle<vtkm::Id> incomingDependentVolume;
rp.dequeue(ingid, incomingDependentVolume);
// TODO/FIXME: Replace with something more efficient
vtkm::Id numSupernodesToProcess =
b->HierarchicalContourTree.FirstSupernodePerIteration[roundNo].ReadPortal().Get(0);
/*
std::cout << "Block " << b->BlockNo << " Round " << roundNo << ": " << numSupernodesToProcess << " entries to process" << std::endl;
vtkm::worklet::contourtree_augmented::PrintIndices(
"Intrinsic Volume (B)", b->IntrinsicVolume, -1, std::cout);
vtkm::worklet::contourtree_augmented::PrintIndices(
"Dependent Volume (B)", b->DependentVolume, -1, std::cout);
*/
auto intrinsicVolumeView =
make_ArrayHandleView(b->IntrinsicVolume, 0, numSupernodesToProcess);
auto incomingIntrinsicVolumeView =
make_ArrayHandleView(incomingIntrinsicVolume, 0, numSupernodesToProcess);
vtkm::cont::ArrayHandle<vtkm::Id> tempSum;
vtkm::cont::Algorithm::Transform(
intrinsicVolumeView, incomingIntrinsicVolumeView, tempSum, vtkm::Sum());
vtkm::cont::Algorithm::Copy(tempSum, intrinsicVolumeView);
auto dependentVolumeView =
make_ArrayHandleView(b->DependentVolume, 0, numSupernodesToProcess);
auto incomingDependentVolumeView =
make_ArrayHandleView(incomingDependentVolume, 0, numSupernodesToProcess);
vtkm::cont::Algorithm::Transform(
dependentVolumeView, incomingDependentVolumeView, tempSum, vtkm::Sum());
vtkm::cont::Algorithm::Copy(tempSum, dependentVolumeView);
/*
vtkm::worklet::contourtree_augmented::PrintIndices(
"Intrinsic Volume (A)", b->IntrinsicVolume, -1, std::cout);
vtkm::worklet::contourtree_augmented::PrintIndices(
"Dependent Volume (A)", b->DependentVolume, -1, std::cout);
*/
}
}
for (int cc = 0; cc < rp.out_link().size(); ++cc)
{
auto target = rp.out_link().target(cc);
if (target.gid != selfid)
{
rp.enqueue(target, b->BlockNo);
rp.enqueue(target, b->IntrinsicVolume);
rp.enqueue(target, b->DependentVolume);
}
}
}
};
void TestHierarchicalHyperSweeper()
{
using vtkm::cont::testing::Testing;
using ContourTreeDataFieldType = vtkm::FloatDefault;
const int numBlocks = 4;
const char* filenames[numBlocks] = { "misc/8x9test_HierarchicalAugmentedTree_Block0.dat",
"misc/8x9test_HierarchicalAugmentedTree_Block1.dat",
"misc/8x9test_HierarchicalAugmentedTree_Block2.dat",
"misc/8x9test_HierarchicalAugmentedTree_Block3.dat" };
vtkm::Id3 globalSize{ 9, 8, 1 };
vtkm::Id3 blocksPerDim{ 2, 2, 1 };
vtkm::Id3 sizes[numBlocks] = { { 5, 4, 1 }, { 5, 5, 1 }, { 5, 4, 1 }, { 5, 5, 1 } };
vtkm::Id3 origins[numBlocks] = { { 0, 0, 0 }, { 0, 3, 0 }, { 4, 0, 0 }, { 4, 3, 0 } };
auto blockIndicesAH = vtkm::cont::make_ArrayHandle(
{ vtkm::Id3{ 0, 0, 0 }, vtkm::Id3{ 0, 1, 0 }, vtkm::Id3{ 1, 0, 0 }, vtkm::Id3{ 1, 1, 0 } });
auto originsAH = vtkm::cont::make_ArrayHandle(
{ vtkm::Id3{ 0, 0, 0 }, vtkm::Id3{ 0, 3, 0 }, vtkm::Id3{ 4, 0, 0 }, vtkm::Id3{ 4, 3, 0 } });
auto sizesAH = vtkm::cont::make_ArrayHandle(
{ vtkm::Id3{ 5, 4, 1 }, vtkm::Id3{ 5, 5, 1 }, vtkm::Id3{ 5, 4, 1 }, vtkm::Id3{ 5, 5, 1 } });
vtkm::worklet::contourtree_distributed::SpatialDecomposition spatialDecomp(
blocksPerDim, globalSize, blockIndicesAH, originsAH, sizesAH);
// Load trees
vtkm::worklet::contourtree_distributed::HierarchicalContourTree<vtkm::FloatDefault>
hct[numBlocks];
for (vtkm::Id blockNo = 0; blockNo < numBlocks; ++blockNo)
{
hct[blockNo].Load(Testing::DataPath(filenames[blockNo]).c_str());
std::cout << hct[blockNo].DebugPrint("AfterLoad", __FILE__, __LINE__);
}
// Create and add DIY blocks
auto comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
vtkm::Id rank = comm.rank();
vtkmdiy::Master master(comm,
1, // Use 1 thread, VTK-M will do the treading
-1 // All block in memory
);
// Set up connectivity
using RegularDecomposer = vtkmdiy::RegularDecomposer<vtkmdiy::DiscreteBounds>;
RegularDecomposer::BoolVector shareFace(3, true);
RegularDecomposer::BoolVector wrap(3, false);
RegularDecomposer::CoordinateVector ghosts(3, 1);
RegularDecomposer::DivisionsVector diyDivisions{ 2, 2, 1 }; // HARDCODED FOR TEST
int numDims = static_cast<int>(globalSize[2] > 1 ? 3 : 2);
RegularDecomposer decomposer(numDims,
spatialDecomp.GetVTKmDIYBounds(),
static_cast<int>(spatialDecomp.GetGlobalNumberOfBlocks()),
shareFace,
wrap,
ghosts,
diyDivisions);
// ... coordinates of local blocks
auto localBlockIndicesPortal = spatialDecomp.LocalBlockIndices.ReadPortal();
std::vector<int> vtkmdiyLocalBlockGids(numBlocks);
for (vtkm::Id bi = 0; bi < numBlocks; bi++)
{
RegularDecomposer::DivisionsVector diyCoords(static_cast<size_t>(numDims));
auto currentCoords = localBlockIndicesPortal.Get(bi);
for (vtkm::IdComponent d = 0; d < numDims; ++d)
{
diyCoords[d] = static_cast<int>(currentCoords[d]);
}
vtkmdiyLocalBlockGids[static_cast<size_t>(bi)] =
RegularDecomposer::coords_to_gid(diyCoords, diyDivisions);
}
// Define which blocks live on which rank so that vtkmdiy can manage them
vtkmdiy::DynamicAssigner assigner(
comm, comm.size(), static_cast<int>(spatialDecomp.GetGlobalNumberOfBlocks()));
for (vtkm::Id bi = 0; bi < numBlocks; bi++)
{
assigner.set_rank(static_cast<int>(rank),
static_cast<int>(vtkmdiyLocalBlockGids[static_cast<size_t>(bi)]));
}
vtkmdiy::fix_links(master, assigner);
HyperSweepBlock<ContourTreeDataFieldType>* localHyperSweeperBlocks[numBlocks];
for (vtkm::Id blockNo = 0; blockNo < numBlocks; ++blockNo)
{
localHyperSweeperBlocks[blockNo] = new HyperSweepBlock<ContourTreeDataFieldType>(
blockNo, origins[blockNo], sizes[blockNo], globalSize, hct[blockNo]);
std::cout << blockNo << " -> " << vtkmdiyLocalBlockGids[blockNo] << std::endl;
master.add(
vtkmdiyLocalBlockGids[blockNo], localHyperSweeperBlocks[blockNo], new vtkmdiy::Link());
}
master.foreach ([](HyperSweepBlock<ContourTreeDataFieldType>* b,
const vtkmdiy::Master::ProxyWithLink&) {
// Create HyperSweeper
std::cout << "Block " << b->BlockNo << std::endl;
std::cout << b->HierarchicalContourTree.DebugPrint(
"Before initializing HyperSweeper", __FILE__, __LINE__);
vtkm::worklet::contourtree_distributed::HierarchicalHyperSweeper<vtkm::Id,
ContourTreeDataFieldType>
hyperSweeper(b->BlockNo, b->HierarchicalContourTree, b->IntrinsicVolume, b->DependentVolume);
std::cout << "Block " << b->BlockNo << std::endl;
std::cout << b->HierarchicalContourTree.DebugPrint(
"After initializing HyperSweeper", __FILE__, __LINE__);
// Create mesh and initialize vertex counts
vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler idRelabeler{ b->Origin,
b->Size,
b->GlobalSize };
if (b->GlobalSize[2] <= 1)
{
vtkm::worklet::contourtree_augmented::DataSetMeshTriangulation2DFreudenthal mesh(
vtkm::Id2{ b->Size[0], b->Size[1] });
hyperSweeper.InitializeIntrinsicVertexCount(
b->HierarchicalContourTree, mesh, idRelabeler, b->IntrinsicVolume);
}
else
{
// TODO/FIXME: For getting owned vertices, it should not make a difference if marching
// cubes or not. Verify.
vtkm::worklet::contourtree_augmented::DataSetMeshTriangulation3DFreudenthal mesh(b->Size);
hyperSweeper.InitializeIntrinsicVertexCount(
b->HierarchicalContourTree, mesh, idRelabeler, b->IntrinsicVolume);
}
std::cout << "Block " << b->BlockNo << std::endl;
std::cout << b->HierarchicalContourTree.DebugPrint(
"After initializing intrinsic vertex count", __FILE__, __LINE__);
// Initialize dependentVolume by copy from intrinsicVolume
vtkm::cont::Algorithm::Copy(b->IntrinsicVolume, b->DependentVolume);
// Perform the local hypersweep
hyperSweeper.LocalHyperSweep();
std::cout << "Block " << b->BlockNo << std::endl;
std::cout << b->HierarchicalContourTree.DebugPrint(
"After local hypersweep", __FILE__, __LINE__);
});
// Reduce
// partners for merge over regular block grid
vtkmdiy::RegularSwapPartners partners(
decomposer, // domain decomposition
2, // radix of k-ary reduction.
true // contiguous: true=distance doubling, false=distance halving
);
vtkmdiy::reduce(
master, assigner, partners, CobmineHyperSweepBlockFunctor<ContourTreeDataFieldType>{});
// Print
vtkm::Id totalVolume = globalSize[0] * globalSize[1] * globalSize[2];
master.foreach ([&totalVolume](HyperSweepBlock<ContourTreeDataFieldType>* b,
const vtkmdiy::Master::ProxyWithLink&) {
std::cout << "Block " << b->BlockNo << std::endl;
std::cout << "=========" << std::endl;
vtkm::worklet::contourtree_augmented::PrintHeader(b->IntrinsicVolume.GetNumberOfValues(),
std::cout);
vtkm::worklet::contourtree_augmented::PrintIndices(
"Intrinsic Volume", b->IntrinsicVolume, -1, std::cout);
vtkm::worklet::contourtree_augmented::PrintIndices(
"Dependent Volume", b->DependentVolume, -1, std::cout);
std::cout << b->HierarchicalContourTree.DebugPrint(
"Called from DumpVolumes", __FILE__, __LINE__);
std::cout << b->HierarchicalContourTree.DumpVolumes(
totalVolume, b->IntrinsicVolume, b->DependentVolume);
});
// Clean-up
for (auto b : localHyperSweeperBlocks)
{
delete b;
}
}
void TestContourTreeUniformDistributed()
{
/*
using vtkm::cont::testing::Testing;
TestContourTreeMeshCombine<vtkm::FloatDefault>(
Testing::DataPath("misc/5x6_7_MC_Rank0_Block0_Round1_BeforeCombineMesh1.ctm"),
Testing::DataPath("misc/5x6_7_MC_Rank0_Block0_Round1_BeforeCombineMesh2.ctm"),
Testing::RegressionImagePath("5x6_7_MC_Rank0_Block0_Round1_CombinedMesh.ctm"));
*/
TestHierarchicalHyperSweeper();
}
} // anonymous namespace