Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-08 13:23:51 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Merge topic 'parallel_pathline2'

Browse Source 
fcd0986c5 Merge branch 'parallel_pathline2' of https://gitlab.kitware.com/dpugmire/vtk-m into parallel_pathline2
ed7358ed2 Kick the dashboard.
6367e4680 Bug fix when pathline rank has 0 ds.
6608a38ee Debugging the dashboard
b4768f06b Bug fix. Must copy particles.
4773b51ca Debug test fail...
982b693ab Forgot to add hxx file.
aa713b565 Support for dist-memory pathlines.

Acked-by: Kitware Robot <kwrobot@kitware.com>
Merge-request: !2399
This commit is contained in:
Dave Pugmire 2021-02-06 02:56:35 +00:00 committed by Kitware Robot
commit 037bd24e7f
17 changed files with 801 additions and 342 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
vtkm/cont/ParticleArrayCopy.h
View File

@ -31,6 +31,17 @@ VTKM_ALWAYS_EXPORT inline void ParticleArrayCopy(
vtkm::cont::ArrayHandle<vtkm::Vec3f, vtkm::cont::StorageTagBasic>& outPos,
bool CopyTerminatedOnly = false);
/// \brief Copy fields in vtkm::Particle to standard types.
///
/// Given a std::vector of \c ArrayHandle of vtkm::Particle, this function copies the
/// position field into an \c ArrayHandle of \c Vec3f objects.
///
template <typename ParticleType>
VTKM_ALWAYS_EXPORT inline void ParticleArrayCopy(
const std::vector<vtkm::cont::ArrayHandle<ParticleType, vtkm::cont::StorageTagBasic>>& inputs,
vtkm::cont::ArrayHandle<vtkm::Vec3f, vtkm::cont::StorageTagBasic>& outPos);
/// \brief Copy all fields in vtkm::Particle to standard types.
///
/// Given an \c ArrayHandle of vtkm::Particle, this function copies the

23
vtkm/cont/ParticleArrayCopy.hxx
View File

@ -87,6 +87,29 @@ VTKM_ALWAYS_EXPORT inline void ParticleArrayCopy(
vtkm::cont::ArrayCopy(posTrn, outPos);
}
template <typename ParticleType>
VTKM_ALWAYS_EXPORT inline void ParticleArrayCopy(
const std::vector<vtkm::cont::ArrayHandle<ParticleType, vtkm::cont::StorageTagBasic>>& inputs,
vtkm::cont::ArrayHandle<vtkm::Vec3f, vtkm::cont::StorageTagBasic>& outPos)
{
vtkm::Id numParticles = 0;
for (const auto& v : inputs)
numParticles += v.GetNumberOfValues();
outPos.Allocate(numParticles);
vtkm::Id idx = 0;
for (const auto& v : inputs)
{
auto posTrn =
vtkm::cont::make_ArrayHandleTransform(v, detail::ExtractPositionFunctor<ParticleType>());
vtkm::Id n = posTrn.GetNumberOfValues();
vtkm::cont::Algorithm::CopySubRange(posTrn, 0, n, outPos, idx);
idx += n;
}
}
/// \brief Copy all fields in vtkm::Particle to standard types.
///
/// Given an \c ArrayHandle of vtkm::Particle, this function copies the

44
vtkm/cont/testing/UnitTestParticleArrayCopy.cxx
View File

@ -29,7 +29,7 @@ void TestParticleArrayCopy()
particles.push_back(vtkm::Particle(vtkm::Vec3f(x, y, z), i));
}
for (int i = 0; i < 2; i++)
for (vtkm::Id i = 0; i < 2; i++)
{
auto particleAH = vtkm::cont::make_ArrayHandle(particles, vtkm::CopyFlag::Off);
@ -64,11 +64,51 @@ void TestParticleArrayCopy()
VTKM_TEST_ASSERT(p.Pos == pt, "Positions do not match");
VTKM_TEST_ASSERT(p.ID == ids.ReadPortal().Get(j), "IDs do not match");
VTKM_TEST_ASSERT(p.NumSteps == steps.ReadPortal().Get(j), "Steps do not match");
VTKM_TEST_ASSERT(p.Status == status.ReadPortal().Get(j), "Steps do not match");
VTKM_TEST_ASSERT(p.Status == status.ReadPortal().Get(j), "Status do not match");
VTKM_TEST_ASSERT(p.Time == ptime.ReadPortal().Get(j), "Times do not match");
}
}
}
//Test copying a vector of ArrayHandles.
std::vector<vtkm::cont::ArrayHandle<vtkm::Particle>> particleVec;
vtkm::Id totalNumParticles = 0;
vtkm::Id pid = 0;
for (vtkm::Id i = 0; i < 4; i++)
{
vtkm::Id n = 5 + i;
std::vector<vtkm::Particle> vec;
for (vtkm::Id j = 0; j < n; j++)
{
auto x = dist(generator);
auto y = dist(generator);
auto z = dist(generator);
vec.push_back(vtkm::Particle(vtkm::Vec3f(x, y, z), pid));
pid++;
}
auto ah = vtkm::cont::make_ArrayHandle(vec, vtkm::CopyFlag::On);
particleVec.push_back(ah);
totalNumParticles += ah.GetNumberOfValues();
}
vtkm::cont::ArrayHandle<vtkm::Vec3f> res;
vtkm::cont::ParticleArrayCopy<vtkm::Particle>(particleVec, res);
VTKM_TEST_ASSERT(res.GetNumberOfValues() == totalNumParticles, "Wrong number of particles");
vtkm::Id resIdx = 0;
auto resPortal = res.ReadPortal();
for (const auto& v : particleVec)
{
vtkm::Id n = v.GetNumberOfValues();
auto portal = v.ReadPortal();
for (vtkm::Id i = 0; i < n; i++)
{
auto p = portal.Get(i);
auto pRes = resPortal.Get(resIdx);
VTKM_TEST_ASSERT(p.Pos == pRes, "Positions do not match");
resIdx++;
}
}
}
int UnitTestParticleArrayCopy(int argc, char* argv[])

17
vtkm/filter/ParticleAdvection.hxx
View File

@ -46,28 +46,33 @@ inline VTKM_CONT vtkm::cont::PartitionedDataSet ParticleAdvection::PrepareForExe
const vtkm::cont::PartitionedDataSet& input,
const vtkm::filter::PolicyBase<DerivedPolicy>&)
{
if (this->GetUseCoordinateSystemAsField())
throw vtkm::cont::ErrorFilterExecution("Coordinate system as field not supported");
if (this->Seeds.GetNumberOfValues() == 0)
throw vtkm::cont::ErrorFilterExecution("No seeds provided.");
std::string activeField = this->GetActiveFieldName();
vtkm::filter::particleadvection::BoundsMap boundsMap(input);
using DataSetIntegratorType = vtkm::filter::particleadvection::DataSetIntegrator;
std::vector<DataSetIntegratorType> dsi;
using DSIType = vtkm::filter::particleadvection::DataSetIntegrator;
std::vector<DSIType> dsi;
for (vtkm::Id i = 0; i < input.GetNumberOfPartitions(); i++)
{
vtkm::Id blockId = boundsMap.GetLocalBlockId(i);
dsi.push_back(DataSetIntegratorType(input.GetPartition(i), blockId, activeField));
auto ds = input.GetPartition(i);
if (!ds.HasPointField(activeField))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
dsi.push_back(DSIType(ds, blockId, activeField));
}
using AlgorithmType = vtkm::filter::particleadvection::ParticleAdvectionAlgorithm;
using ThreadedAlgorithmType = vtkm::filter::particleadvection::ParticleAdvectionThreadedAlgorithm;
if (this->UseThreadedAlgorithm)
return vtkm::filter::particleadvection::RunAlgo<ThreadedAlgorithmType>(
if (this->GetUseThreadedAlgorithm())
return vtkm::filter::particleadvection::RunAlgo<DSIType, ThreadedAlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
else
return vtkm::filter::particleadvection::RunAlgo<AlgorithmType>(
return vtkm::filter::particleadvection::RunAlgo<DSIType, AlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
}

27
vtkm/filter/Pathline.h
View File

@ -38,7 +38,14 @@ public:
void SetNextTime(vtkm::FloatDefault t) { this->NextTime = t; }
VTKM_CONT
void SetNextDataSet(const vtkm::cont::DataSet& ds) { this->NextDataSet = ds; }
void SetNextDataSet(const vtkm::cont::DataSet& ds)
{
this->NextDataSet = vtkm::cont::PartitionedDataSet(ds);
}
VTKM_CONT
void SetNextDataSet(const vtkm::cont::PartitionedDataSet& pds) { this->NextDataSet = pds; }
VTKM_CONT
void SetStepSize(vtkm::FloatDefault s) { this->StepSize = s; }
@ -49,11 +56,15 @@ public:
VTKM_CONT
void SetSeeds(vtkm::cont::ArrayHandle<vtkm::Particle>& seeds);
template <typename T, typename StorageType, typename DerivedPolicy>
VTKM_CONT vtkm::cont::DataSet DoExecute(
const vtkm::cont::DataSet& input,
const vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, StorageType>& field,
const vtkm::filter::FieldMetadata& fieldMeta,
VTKM_CONT
bool GetUseThreadedAlgorithm() { return this->UseThreadedAlgorithm; }
VTKM_CONT
void SetUseThreadedAlgorithm(bool val) { this->UseThreadedAlgorithm = val; }
template <typename DerivedPolicy>
vtkm::cont::PartitionedDataSet PrepareForExecution(
const vtkm::cont::PartitionedDataSet& input,
const vtkm::filter::PolicyBase<DerivedPolicy>& policy);
template <typename DerivedPolicy>
@ -62,13 +73,13 @@ public:
vtkm::filter::PolicyBase<DerivedPolicy> policy);
private:
vtkm::worklet::Streamline Worklet;
vtkm::FloatDefault StepSize;
vtkm::FloatDefault PreviousTime;
vtkm::FloatDefault NextTime;
vtkm::cont::DataSet NextDataSet;
vtkm::cont::PartitionedDataSet NextDataSet;
vtkm::Id NumberOfSteps;
vtkm::cont::ArrayHandle<vtkm::Particle> Seeds;
bool UseThreadedAlgorithm;
};
}
} // namespace vtkm::filter

84
vtkm/filter/Pathline.hxx
View File

@ -11,15 +11,14 @@
#ifndef vtk_m_filter_Pathline_hxx
#define vtk_m_filter_Pathline_hxx
#include <vtkm/filter/Pathline.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ArrayHandleIndex.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/worklet/particleadvection/Field.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
#include <vtkm/worklet/particleadvection/Particles.h>
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
#include <vtkm/worklet/particleadvection/TemporalGridEvaluators.h>
#include <vtkm/filter/particleadvection/BoundsMap.h>
#include <vtkm/filter/particleadvection/DataSetIntegrator.h>
#include <vtkm/filter/particleadvection/StreamlineAlgorithm.h>
namespace vtkm
{
@ -29,7 +28,7 @@ namespace filter
//-----------------------------------------------------------------------------
inline VTKM_CONT Pathline::Pathline()
: vtkm::filter::FilterDataSetWithField<Pathline>()
, Worklet()
, UseThreadedAlgorithm(false)
{
}
@ -40,58 +39,45 @@ inline VTKM_CONT void Pathline::SetSeeds(vtkm::cont::ArrayHandle<vtkm::Particle>
}
//-----------------------------------------------------------------------------
template <typename T, typename StorageType, typename DerivedPolicy>
inline VTKM_CONT vtkm::cont::DataSet Pathline::DoExecute(
const vtkm::cont::DataSet& input,
const vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, StorageType>& field,
const vtkm::filter::FieldMetadata& fieldMeta,
template <typename DerivedPolicy>
inline VTKM_CONT vtkm::cont::PartitionedDataSet Pathline::PrepareForExecution(
const vtkm::cont::PartitionedDataSet& input,
const vtkm::filter::PolicyBase<DerivedPolicy>&)
{
//Check for some basics.
if (this->GetUseCoordinateSystemAsField())
throw vtkm::cont::ErrorFilterExecution("Coordinate system as field not supported");
if (this->Seeds.GetNumberOfValues() == 0)
{
throw vtkm::cont::ErrorFilterExecution("No seeds provided.");
}
if (this->NextDataSet.GetNumberOfPartitions() != input.GetNumberOfPartitions())
throw vtkm::cont::ErrorFilterExecution("Number of partitions do not match");
if (!(this->PreviousTime < this->NextTime))
throw vtkm::cont::ErrorFilterExecution("Previous time must be less than Next time.");
const vtkm::cont::DynamicCellSet& cells = input.GetCellSet();
const vtkm::cont::DynamicCellSet& cells2 = this->NextDataSet.GetCellSet();
const vtkm::cont::CoordinateSystem& coords =
input.GetCoordinateSystem(this->GetActiveCoordinateSystemIndex());
const vtkm::cont::CoordinateSystem& coords2 =
this->NextDataSet.GetCoordinateSystem(this->GetActiveCoordinateSystemIndex());
std::string activeField = this->GetActiveFieldName();
vtkm::filter::particleadvection::BoundsMap boundsMap(input);
using DSIType = vtkm::filter::particleadvection::TemporalDataSetIntegrator;
std::vector<DSIType> dsi;
auto field2 = vtkm::cont::Cast<vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, StorageType>>(
this->NextDataSet.GetField(this->GetActiveFieldName()).GetData());
if (!fieldMeta.IsPointField())
for (vtkm::Id i = 0; i < input.GetNumberOfPartitions(); i++)
{
throw vtkm::cont::ErrorFilterExecution("Point field expected.");
vtkm::Id blockId = boundsMap.GetLocalBlockId(i);
auto dsPrev = input.GetPartition(i);
auto dsNext = this->NextDataSet.GetPartition(i);
if (!dsPrev.HasPointField(activeField) || !dsNext.HasPointField(activeField))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
dsi.push_back(
DSIType(dsPrev, this->PreviousTime, dsNext, this->NextTime, blockId, activeField));
}
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, StorageType>;
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::particleadvection::TemporalGridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
using AlgorithmType = vtkm::filter::particleadvection::PathlineAlgorithm;
using ThreadedAlgorithmType = vtkm::filter::particleadvection::PathlineThreadedAlgorithm;
FieldType velocities(field);
FieldType velocities2(field2);
GridEvalType eval(
coords, cells, velocities, this->PreviousTime, coords2, cells2, velocities2, this->NextTime);
RK4Type rk4(eval, this->StepSize);
vtkm::worklet::Streamline streamline;
vtkm::worklet::StreamlineResult<vtkm::Particle> res;
vtkm::cont::ArrayHandle<vtkm::Particle> seedArray;
vtkm::cont::ArrayCopy(this->Seeds, seedArray);
res = Worklet.Run(rk4, seedArray, this->NumberOfSteps);
vtkm::cont::DataSet outData;
vtkm::cont::CoordinateSystem outputCoords("coordinates", res.Positions);
outData.SetCellSet(res.PolyLines);
outData.AddCoordinateSystem(outputCoords);
return outData;
if (this->GetUseThreadedAlgorithm())
return vtkm::filter::particleadvection::RunAlgo<DSIType, ThreadedAlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
else
return vtkm::filter::particleadvection::RunAlgo<DSIType, AlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
}
//-----------------------------------------------------------------------------

17
vtkm/filter/Streamline.hxx
View File

@ -45,28 +45,33 @@ inline VTKM_CONT vtkm::cont::PartitionedDataSet Streamline::PrepareForExecution(
const vtkm::cont::PartitionedDataSet& input,
const vtkm::filter::PolicyBase<DerivedPolicy>&)
{
if (this->GetUseCoordinateSystemAsField())
throw vtkm::cont::ErrorFilterExecution("Coordinate system as field not supported");
if (this->Seeds.GetNumberOfValues() == 0)
throw vtkm::cont::ErrorFilterExecution("No seeds provided.");
std::string activeField = this->GetActiveFieldName();
vtkm::filter::particleadvection::BoundsMap boundsMap(input);
using DataSetIntegratorType = vtkm::filter::particleadvection::DataSetIntegrator;
std::vector<DataSetIntegratorType> dsi;
using DSIType = vtkm::filter::particleadvection::DataSetIntegrator;
std::vector<DSIType> dsi;
for (vtkm::Id i = 0; i < input.GetNumberOfPartitions(); i++)
{
vtkm::Id blockId = boundsMap.GetLocalBlockId(i);
dsi.push_back(DataSetIntegratorType(input.GetPartition(i), blockId, activeField));
auto ds = input.GetPartition(i);
if (!ds.HasPointField(activeField))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
dsi.push_back(DSIType(ds, blockId, activeField));
}
using AlgorithmType = vtkm::filter::particleadvection::StreamlineAlgorithm;
using ThreadedAlgorithmType = vtkm::filter::particleadvection::StreamlineThreadedAlgorithm;
if (this->UseThreadedAlgorithm)
return vtkm::filter::particleadvection::RunAlgo<ThreadedAlgorithmType>(
if (this->GetUseThreadedAlgorithm())
return vtkm::filter::particleadvection::RunAlgo<DSIType, ThreadedAlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
else
return vtkm::filter::particleadvection::RunAlgo<AlgorithmType>(
return vtkm::filter::particleadvection::RunAlgo<DSIType, AlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
}

396
vtkm/filter/particleadvection/AdvectorBaseAlgorithm.h
View File

@ -11,7 +11,10 @@
#ifndef vtk_m_filter_particleadvection_AdvectorBaseAlgorithm_h
#define vtk_m_filter_particleadvection_AdvectorBaseAlgorithm_h
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayHandlePermutation.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/cont/ParticleArrayCopy.h>
#include <vtkm/filter/particleadvection/DataSetIntegrator.h>
#include <vtkm/filter/particleadvection/ParticleMessenger.h>
@ -24,9 +27,183 @@ namespace filter
{
namespace particleadvection
{
using DataSetIntegratorType = vtkm::filter::particleadvection::DataSetIntegrator;
template <typename AlgorithmType>
namespace internal
{
//Helper class to store the different result types.
template <typename ResultType>
class ResultHelper;
//Specialization for ParticleAdvectionResult
using PAType = vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>;
template <>
class ResultHelper<PAType>
{
public:
static void Store(std::map<vtkm::Id, std::vector<PAType>>& results,
const PAType& res,
vtkm::Id blockId,
const std::vector<vtkm::Id>& indices)
{
if (indices.empty())
return;
//Selected out the terminated particles and store them.
vtkm::cont::ArrayHandle<vtkm::Particle> termParticles;
auto indicesAH = vtkm::cont::make_ArrayHandle(indices, vtkm::CopyFlag::Off);
auto perm = vtkm::cont::make_ArrayHandlePermutation(indicesAH, res.Particles);
vtkm::cont::Algorithm::Copy(perm, termParticles);
PAType termRes(termParticles);
results[blockId].push_back(termRes);
}
static vtkm::cont::PartitionedDataSet GetOutput(
const std::map<vtkm::Id, std::vector<PAType>>& results)
{
vtkm::cont::PartitionedDataSet output;
for (const auto& it : results)
{
std::size_t nResults = it.second.size();
if (nResults == 0)
continue;
std::vector<vtkm::cont::ArrayHandle<vtkm::Particle>> allParticles;
allParticles.reserve(static_cast<std::size_t>(nResults));
for (const auto& res : it.second)
allParticles.push_back(res.Particles);
vtkm::cont::ArrayHandle<vtkm::Vec3f> pts;
vtkm::cont::ParticleArrayCopy(allParticles, pts);
vtkm::cont::DataSet ds;
vtkm::cont::CoordinateSystem outCoords("coordinates", pts);
ds.AddCoordinateSystem(outCoords);
//Create vertex cell set
vtkm::Id numPoints = pts.GetNumberOfValues();
vtkm::cont::CellSetSingleType<> cells;
vtkm::cont::ArrayHandleIndex conn(numPoints);
vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
vtkm::cont::ArrayCopy(conn, connectivity);
cells.Fill(numPoints, vtkm::CELL_SHAPE_VERTEX, 1, connectivity);
ds.SetCellSet(cells);
output.AppendPartition(ds);
}
return output;
}
};
//Specialization for StreamlineResult
using SLType = vtkm::worklet::StreamlineResult<vtkm::Particle>;
template <>
class ResultHelper<SLType>
{
public:
static void Store(std::map<vtkm::Id, std::vector<SLType>>& results,
const SLType& res,
vtkm::Id blockId,
const std::vector<vtkm::Id>& vtkmNotUsed(indices))
{
results[blockId].push_back(res);
}
static vtkm::cont::PartitionedDataSet GetOutput(
const std::map<vtkm::Id, std::vector<SLType>>& results)
{
vtkm::cont::PartitionedDataSet output;
for (const auto& it : results)
{
std::size_t nResults = it.second.size();
if (nResults == 0)
continue;
vtkm::cont::DataSet ds;
//Easy case with one result.
if (nResults == 1)
{
ds.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", it.second[0].Positions));
ds.SetCellSet(it.second[0].PolyLines);
}
else
{
//Append all the results into one data set.
vtkm::cont::ArrayHandle<vtkm::Vec3f> appendPts;
std::vector<vtkm::Id> posOffsets(nResults);
const auto& res0 = it.second[0];
vtkm::Id totalNumCells = res0.PolyLines.GetNumberOfCells();
vtkm::Id totalNumPts = res0.Positions.GetNumberOfValues();
posOffsets[0] = 0;
for (std::size_t i = 1; i < nResults; i++)
{
const auto& res = it.second[i];
posOffsets[i] = totalNumPts;
totalNumPts += res.Positions.GetNumberOfValues();
totalNumCells += res.PolyLines.GetNumberOfCells();
}
//Append all the points together.
appendPts.Allocate(totalNumPts);
for (std::size_t i = 0; i < nResults; i++)
{
const auto& res = it.second[i];
// copy all values into appendPts starting at offset.
vtkm::cont::Algorithm::CopySubRange(
res.Positions, 0, res.Positions.GetNumberOfValues(), appendPts, posOffsets[i]);
}
vtkm::cont::CoordinateSystem outputCoords =
vtkm::cont::CoordinateSystem("coordinates", appendPts);
ds.AddCoordinateSystem(outputCoords);
//Create polylines.
std::vector<vtkm::Id> numPtsPerCell(static_cast<std::size_t>(totalNumCells));
std::size_t off = 0;
for (std::size_t i = 0; i < nResults; i++)
{
const auto& res = it.second[i];
vtkm::Id nCells = res.PolyLines.GetNumberOfCells();
for (vtkm::Id j = 0; j < nCells; j++)
numPtsPerCell[off++] = static_cast<vtkm::Id>(res.PolyLines.GetNumberOfPointsInCell(j));
}
auto numPointsPerCellArray =
vtkm::cont::make_ArrayHandle(numPtsPerCell, vtkm::CopyFlag::Off);
vtkm::cont::ArrayHandle<vtkm::Id> cellIndex;
vtkm::Id connectivityLen =
vtkm::cont::Algorithm::ScanExclusive(numPointsPerCellArray, cellIndex);
vtkm::cont::ArrayHandleIndex connCount(connectivityLen);
vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
vtkm::cont::ArrayCopy(connCount, connectivity);
vtkm::cont::ArrayHandle<vtkm::UInt8> cellTypes;
auto polyLineShape = vtkm::cont::make_ArrayHandleConstant<vtkm::UInt8>(
vtkm::CELL_SHAPE_POLY_LINE, totalNumCells);
vtkm::cont::ArrayCopy(polyLineShape, cellTypes);
auto offsets = vtkm::cont::ConvertNumIndicesToOffsets(numPointsPerCellArray);
vtkm::cont::CellSetExplicit<> polyLines;
polyLines.Fill(totalNumPts, cellTypes, connectivity, offsets);
ds.SetCellSet(polyLines);
}
output.AppendPartition(ds);
}
return output;
}
};
};
template <typename DataSetIntegratorType, typename AlgorithmType>
vtkm::cont::PartitionedDataSet RunAlgo(const vtkm::filter::particleadvection::BoundsMap& boundsMap,
const std::vector<DataSetIntegratorType>& dsi,
vtkm::Id numSteps,
@ -45,11 +222,9 @@ vtkm::cont::PartitionedDataSet RunAlgo(const vtkm::filter::particleadvection::Bo
//
// Base class for particle advector
//
template <typename ResultType>
template <typename DataSetIntegratorType, typename ResultType>
class VTKM_ALWAYS_EXPORT AdvectorBaseAlgorithm
{
using DataSetIntegratorType = vtkm::filter::particleadvection::DataSetIntegrator;
public:
AdvectorBaseAlgorithm(const vtkm::filter::particleadvection::BoundsMap& bm,
const std::vector<DataSetIntegratorType>& blocks)
@ -107,7 +282,6 @@ public:
if (GetActiveParticles(v, blockId))
{
const auto& block = this->GetDataSet(blockId);
ResultType r;
block.Advect(v, this->StepSize, this->NumberOfSteps, r);
numTerm = this->UpdateResult(r, blockId);
@ -122,14 +296,16 @@ public:
}
}
inline vtkm::cont::PartitionedDataSet GetOutput();
vtkm::cont::PartitionedDataSet GetOutput() const
{
return internal::ResultHelper<ResultType>::GetOutput(this->Results);
}
protected:
virtual void ClearParticles()
{
this->Active.clear();
this->Inactive.clear();
this->Terminated.clear();
this->ParticleBlockIDsMap.clear();
}
@ -151,69 +327,6 @@ protected:
throw vtkm::cont::ErrorFilterExecution("Bad block");
}
void UpdateResultParticle(vtkm::Particle& p,
std::vector<vtkm::Particle>& I,
std::vector<vtkm::Particle>& T,
std::vector<vtkm::Particle>& A,
std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& idsMapI,
std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& idsMapA) const
{
if (p.Status.CheckTerminate())
T.push_back(p);
else
{
const auto& it = this->ParticleBlockIDsMap.find(p.ID);
VTKM_ASSERT(it != this->ParticleBlockIDsMap.end());
auto currBIDs = it->second;
VTKM_ASSERT(!currBIDs.empty());
std::vector<vtkm::Id> newIDs;
if (p.Status.CheckSpatialBounds() && !p.Status.CheckTookAnySteps())
newIDs.assign(std::next(currBIDs.begin(), 1), currBIDs.end());
else
newIDs = this->BoundsMap.FindBlocks(p.Pos, currBIDs);
//reset the particle status.
p.Status = vtkm::ParticleStatus();
if (newIDs.empty()) //No blocks, we're done.
{
p.Status.SetTerminate();
T.push_back(p);
}
else
{
//If we have more than blockId, we want to minimize communication
//and put any blocks owned by this rank first.
if (newIDs.size() > 1)
{
for (auto idit = newIDs.begin(); idit != newIDs.end(); idit++)
{
vtkm::Id bid = *idit;
if (this->BoundsMap.FindRank(bid) == this->Rank)
{
newIDs.erase(idit);
newIDs.insert(newIDs.begin(), bid);
break;
}
}
}
int dstRank = this->BoundsMap.FindRank(newIDs[0]);
if (dstRank == this->Rank)
{
A.push_back(p);
idsMapA[p.ID] = newIDs;
}
else
{
I.push_back(p);
idsMapI[p.ID] = newIDs;
}
}
}
}
virtual void SetSeedArray(const std::vector<vtkm::Particle>& particles,
const std::vector<std::vector<vtkm::Id>>& blockIds)
{
@ -277,58 +390,133 @@ protected:
virtual void UpdateActive(const std::vector<vtkm::Particle>& particles,
const std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& idsMap)
{
VTKM_ASSERT(particles.size() == idsMap.size());
if (particles.empty())
return;
this->Active.insert(this->Active.end(), particles.begin(), particles.end());
for (const auto& it : idsMap)
this->ParticleBlockIDsMap[it.first] = it.second;
this->Update(this->Active, particles, idsMap);
}
virtual void UpdateInactive(const std::vector<vtkm::Particle>& particles,
const std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& idsMap)
{
VTKM_ASSERT(particles.size() == idsMap.size());
if (particles.empty())
return;
this->Update(this->Inactive, particles, idsMap);
}
this->Inactive.insert(this->Inactive.end(), particles.begin(), particles.end());
void Update(std::vector<vtkm::Particle>& arr,
const std::vector<vtkm::Particle>& particles,
const std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& idsMap)
{
VTKM_ASSERT(particles.size() == idsMap.size());
arr.insert(arr.end(), particles.begin(), particles.end());
for (const auto& it : idsMap)
this->ParticleBlockIDsMap[it.first] = it.second;
}
virtual void UpdateTerminated(const std::vector<vtkm::Particle>& particles, vtkm::Id blockId)
void UpdateTerminated(const vtkm::cont::ArrayHandle<vtkm::Particle>& particles,
const std::vector<vtkm::Id>& idxTerm)
{
if (particles.empty())
return;
for (const auto& t : particles)
this->ParticleBlockIDsMap.erase(t.ID);
auto& it = this->Terminated[blockId];
it.insert(it.end(), particles.begin(), particles.end());
auto portal = particles.ReadPortal();
for (const auto& idx : idxTerm)
this->ParticleBlockIDsMap.erase(portal.Get(idx).ID);
}
vtkm::Id UpdateResult(const ResultType& res, vtkm::Id blockId)
void ClassifyParticles(const vtkm::cont::ArrayHandle<vtkm::Particle>& particles,
std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& idsMapA,
std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& idsMapI,
std::vector<vtkm::Particle>& A,
std::vector<vtkm::Particle>& I,
std::vector<vtkm::Id>& termIdx) const
{
vtkm::Id n = res.Particles.GetNumberOfValues();
auto portal = res.Particles.ReadPortal();
A.clear();
I.clear();
termIdx.clear();
idsMapI.clear();
idsMapA.clear();
std::unordered_map<vtkm::Id, std::vector<vtkm::Id>> idsMapI, idsMapA;
std::vector<vtkm::Particle> I, T, A;
auto portal = particles.WritePortal();
vtkm::Id n = portal.GetNumberOfValues();
for (vtkm::Id i = 0; i < n; i++)
{
vtkm::Particle p = portal.Get(i);
this->UpdateResultParticle(p, I, T, A, idsMapI, idsMapA);
auto p = portal.Get(i);
if (p.Status.CheckTerminate())
termIdx.push_back(i);
else
{
const auto& it = this->ParticleBlockIDsMap.find(p.ID);
VTKM_ASSERT(it != this->ParticleBlockIDsMap.end());
auto currBIDs = it->second;
VTKM_ASSERT(!currBIDs.empty());
std::vector<vtkm::Id> newIDs;
if (p.Status.CheckSpatialBounds() && !p.Status.CheckTookAnySteps())
newIDs.assign(std::next(currBIDs.begin(), 1), currBIDs.end());
else
newIDs = this->BoundsMap.FindBlocks(p.Pos, currBIDs);
//reset the particle status.
p.Status = vtkm::ParticleStatus();
if (newIDs.empty()) //No blocks, we're done.
{
p.Status.SetTerminate();
termIdx.push_back(i);
}
else
{
//If we have more than blockId, we want to minimize communication
//and put any blocks owned by this rank first.
if (newIDs.size() > 1)
{
for (auto idit = newIDs.begin(); idit != newIDs.end(); idit++)
{
vtkm::Id bid = *idit;
if (this->BoundsMap.FindRank(bid) == this->Rank)
{
newIDs.erase(idit);
newIDs.insert(newIDs.begin(), bid);
break;
}
}
}
int dstRank = this->BoundsMap.FindRank(newIDs[0]);
if (dstRank == this->Rank)
{
A.push_back(p);
idsMapA[p.ID] = newIDs;
}
else
{
I.push_back(p);
idsMapI[p.ID] = newIDs;
}
}
portal.Set(i, p);
}
}
vtkm::Id numTerm = static_cast<vtkm::Id>(T.size());
//Make sure we didn't miss anything. Every particle goes into a single bucket.
VTKM_ASSERT(static_cast<std::size_t>(n) == (A.size() + I.size() + termIdx.size()));
}
vtkm::Id UpdateResult(ResultType& res, vtkm::Id blockId)
{
std::unordered_map<vtkm::Id, std::vector<vtkm::Id>> idsMapI, idsMapA;
std::vector<vtkm::Particle> A, I;
std::vector<vtkm::Id> termIdx;
this->ClassifyParticles(res.Particles, idsMapA, idsMapI, A, I, termIdx);
//Update active, inactive and terminated
this->UpdateActive(A, idsMapA);
this->UpdateInactive(I, idsMapI);
this->UpdateTerminated(T, blockId);
this->StoreResult(res, blockId);
vtkm::Id numTerm = static_cast<vtkm::Id>(termIdx.size());
if (numTerm > 0)
this->UpdateTerminated(res.Particles, termIdx);
internal::ResultHelper<ResultType>::Store(this->Results, res, blockId, termIdx);
return numTerm;
}
@ -348,8 +536,6 @@ protected:
return false;
}
inline void StoreResult(const ResultType& res, vtkm::Id blockId);
//Member data
std::vector<vtkm::Particle> Active;
std::vector<DataSetIntegratorType> Blocks;
@ -364,16 +550,10 @@ protected:
vtkm::FloatDefault StepSize;
vtkm::Id TotalNumParticles;
vtkm::Id TotalNumTerminatedParticles;
std::unordered_map<vtkm::Id, std::vector<vtkm::Particle>> Terminated;
};
}
}
} // namespace vtkm::filter::particleadvection
#ifndef vtk_m_filter_particleadvection_AdvectorBaseAlgorithm_hxx
#include <vtkm/filter/particleadvection/AdvectorBaseAlgorithm.hxx>
#endif
#endif

24
vtkm/filter/particleadvection/AdvectorBaseThreadedAlgorithm.h
View File

@ -22,15 +22,14 @@ namespace filter
namespace particleadvection
{
template <typename ResultType>
class VTKM_ALWAYS_EXPORT AdvectorBaseThreadedAlgorithm : public AdvectorBaseAlgorithm<ResultType>
template <typename DataSetIntegratorType, typename ResultType>
class VTKM_ALWAYS_EXPORT AdvectorBaseThreadedAlgorithm
: public AdvectorBaseAlgorithm<DataSetIntegratorType, ResultType>
{
public:
using DataSetIntegratorType = vtkm::filter::particleadvection::DataSetIntegrator;
AdvectorBaseThreadedAlgorithm(const vtkm::filter::particleadvection::BoundsMap& bm,
const std::vector<DataSetIntegratorType>& blocks)
: AdvectorBaseAlgorithm<ResultType>(bm, blocks)
: AdvectorBaseAlgorithm<DataSetIntegratorType, ResultType>(bm, blocks)
, Done(false)
, WorkerActivate(false)
{
@ -57,7 +56,8 @@ protected:
bool GetActiveParticles(std::vector<vtkm::Particle>& particles, vtkm::Id& blockId) override
{
std::lock_guard<std::mutex> lock(this->Mutex);
bool val = this->AdvectorBaseAlgorithm<ResultType>::GetActiveParticles(particles, blockId);
bool val = this->AdvectorBaseAlgorithm<DataSetIntegratorType, ResultType>::GetActiveParticles(
particles, blockId);
this->WorkerActivate = val;
return val;
}
@ -68,7 +68,8 @@ protected:
if (!particles.empty())
{
std::lock_guard<std::mutex> lock(this->Mutex);
this->AdvectorBaseAlgorithm<ResultType>::UpdateActive(particles, idsMap);
this->AdvectorBaseAlgorithm<DataSetIntegratorType, ResultType>::UpdateActive(particles,
idsMap);
//Let workers know there is new work
this->WorkerActivateCondition.notify_all();
@ -106,7 +107,6 @@ protected:
if (this->GetActiveParticles(v, blockId))
{
const auto& block = this->GetDataSet(blockId);
ResultType r;
block.Advect(v, this->StepSize, this->NumberOfSteps, r);
this->UpdateWorkerResult(blockId, r);
@ -127,11 +127,10 @@ protected:
this->GetWorkerResults(workerResults);
vtkm::Id numTerm = 0;
for (const auto& it : workerResults)
for (auto& it : workerResults)
{
vtkm::Id blockId = it.first;
const auto& results = it.second;
for (const auto& r : results)
for (auto& r : it.second)
numTerm += this->UpdateResult(r, blockId);
}
@ -151,7 +150,8 @@ protected:
{
std::lock_guard<std::mutex> lock(this->Mutex);
return (this->AdvectorBaseAlgorithm<ResultType>::GetBlockAndWait(numLocalTerm) &&
return (this->AdvectorBaseAlgorithm<DataSetIntegratorType, ResultType>::GetBlockAndWait(
numLocalTerm) &&
!this->WorkerActivate && this->WorkerResults.empty());
}

1
vtkm/filter/particleadvection/CMakeLists.txt
View File

@ -14,6 +14,7 @@ set(headers
AdvectorBaseThreadedAlgorithm.h
BoundsMap.h
DataSetIntegrator.h
DataSetIntegrator.hxx
Messenger.h
ParticleAdvectionAlgorithm.h
ParticleMessenger.h

115
vtkm/filter/particleadvection/DataSetIntegrator.h
View File

@ -15,6 +15,7 @@
#include <vtkm/cont/DataSet.h>
#include <vtkm/worklet/ParticleAdvection.h>
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
#include <vtkm/worklet/particleadvection/TemporalGridEvaluators.h>
#include <memory>
#include <vector>
@ -26,31 +27,19 @@ namespace filter
namespace particleadvection
{
class VTKM_ALWAYS_EXPORT DataSetIntegrator
template <typename GridEvalType>
class VTKM_ALWAYS_EXPORT DataSetIntegratorBase
{
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
public:
DataSetIntegrator(const vtkm::cont::DataSet& ds, vtkm::Id id, std::string fieldNm)
: ActiveField(ds.GetField(fieldNm))
, CopySeedArray(false)
DataSetIntegratorBase(bool copySeeds = false, vtkm::Id id = -1)
: CopySeedArray(copySeeds)
, Eval(nullptr)
, ID(id)
{
auto fieldData = this->ActiveField.GetData();
FieldHandle fieldArray;
if (fieldData.IsType<FieldHandle>())
fieldArray = fieldData.AsArrayHandle<FieldHandle>();
else
vtkm::cont::ArrayCopy(fieldData, fieldArray);
this->Eval = std::shared_ptr<GridEvalType>(new GridEvalType(ds, fieldArray));
}
~DataSetIntegratorBase() = default;
vtkm::Id GetID() const { return this->ID; }
void SetCopySeedFlag(bool val) { this->CopySeedArray = val; }
@ -66,46 +55,88 @@ public:
this->DoAdvect(seedArray, rk4, maxSteps, result);
}
private:
protected:
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
using FieldHandleType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
template <typename ResultType>
inline void DoAdvect(vtkm::cont::ArrayHandle<vtkm::Particle>& seeds,
const RK4Type& rk4,
vtkm::Id maxSteps,
ResultType& result) const;
vtkm::cont::Field ActiveField;
FieldHandleType GetFieldHandle(const vtkm::cont::DataSet& ds, const std::string& fieldNm)
{
if (!ds.HasField(fieldNm))
throw vtkm::cont::ErrorFilterExecution("Field " + fieldNm + " not found on dataset.");
FieldHandleType fieldArray;
auto fieldData = ds.GetField(fieldNm).GetData();
if (fieldData.IsType<FieldHandleType>())
fieldArray = fieldData.AsArrayHandle<FieldHandleType>();
else
vtkm::cont::ArrayCopy(
fieldData.ResetTypes<vtkm::TypeListFieldVec3, VTKM_DEFAULT_STORAGE_LIST>(), fieldArray);
return fieldArray;
}
bool CopySeedArray;
std::shared_ptr<GridEvalType> Eval;
vtkm::Id ID;
};
//-----
// Specialization for ParticleAdvection worklet
template <>
inline void DataSetIntegrator::DoAdvect(
vtkm::cont::ArrayHandle<vtkm::Particle>& seeds,
const RK4Type& rk4,
vtkm::Id maxSteps,
vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>& result) const
class VTKM_ALWAYS_EXPORT DataSetIntegrator
: public DataSetIntegratorBase<vtkm::worklet::particleadvection::GridEvaluator<
vtkm::worklet::particleadvection::VelocityField<vtkm::cont::ArrayHandle<vtkm::Vec3f>>>>
{
vtkm::worklet::ParticleAdvection Worklet;
result = Worklet.Run(rk4, seeds, maxSteps);
}
public:
DataSetIntegrator(const vtkm::cont::DataSet& ds, vtkm::Id id, const std::string& fieldNm)
: DataSetIntegratorBase<vtkm::worklet::particleadvection::GridEvaluator<
vtkm::worklet::particleadvection::VelocityField<FieldHandleType>>>(false, id)
{
auto fieldArray = this->GetFieldHandle(ds, fieldNm);
//-----
// Specialization for Streamline worklet
template <>
inline void DataSetIntegrator::DoAdvect(
vtkm::cont::ArrayHandle<vtkm::Particle>& seeds,
const RK4Type& rk4,
vtkm::Id maxSteps,
vtkm::worklet::StreamlineResult<vtkm::Particle>& result) const
using EvalType = vtkm::worklet::particleadvection::GridEvaluator<
vtkm::worklet::particleadvection::VelocityField<FieldHandleType>>;
this->Eval = std::shared_ptr<EvalType>(new EvalType(ds, fieldArray));
}
};
class VTKM_ALWAYS_EXPORT TemporalDataSetIntegrator
: public DataSetIntegratorBase<vtkm::worklet::particleadvection::TemporalGridEvaluator<
vtkm::worklet::particleadvection::VelocityField<vtkm::cont::ArrayHandle<vtkm::Vec3f>>>>
{
vtkm::worklet::Streamline Worklet;
result = Worklet.Run(rk4, seeds, maxSteps);
}
using FieldHandleType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
public:
TemporalDataSetIntegrator(const vtkm::cont::DataSet& ds1,
vtkm::FloatDefault t1,
const vtkm::cont::DataSet& ds2,
vtkm::FloatDefault t2,
vtkm::Id id,
const std::string& fieldNm)
: DataSetIntegratorBase<vtkm::worklet::particleadvection::TemporalGridEvaluator<
vtkm::worklet::particleadvection::VelocityField<FieldHandleType>>>(false, id)
{
auto fieldArray1 = this->GetFieldHandle(ds1, fieldNm);
auto fieldArray2 = this->GetFieldHandle(ds2, fieldNm);
using EvalType = vtkm::worklet::particleadvection::TemporalGridEvaluator<
vtkm::worklet::particleadvection::VelocityField<FieldHandleType>>;
this->Eval =
std::shared_ptr<EvalType>(new EvalType(ds1, t1, fieldArray1, ds2, t2, fieldArray2));
}
};
}
}
} // namespace vtkm::filter::particleadvection
#ifndef vtk_m_filter_particleadvection_DataSetIntegrator_hxx
#include <vtkm/filter/particleadvection/DataSetIntegrator.hxx>
#endif
#endif //vtk_m_filter_DataSetIntegrator_h

74
vtkm/filter/particleadvection/DataSetIntegrator.hxx Normal file
View File

@ -0,0 +1,74 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_DataSetIntegrator_hxx
#define vtk_m_filter_DataSetIntegrator_hxx
namespace vtkm
{
namespace filter
{
namespace particleadvection
{
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<
vtkm::worklet::particleadvection::VelocityField<vtkm::cont::ArrayHandle<vtkm::Vec3f>>>;
using TemporalGridEvalType = vtkm::worklet::particleadvection::TemporalGridEvaluator<
vtkm::worklet::particleadvection::VelocityField<vtkm::cont::ArrayHandle<vtkm::Vec3f>>>;
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
using TemporalRK4Type = vtkm::worklet::particleadvection::RK4Integrator<TemporalGridEvalType>;
//-----
// Specialization for ParticleAdvection worklet
template <>
template <>
inline void DataSetIntegratorBase<GridEvalType>::DoAdvect(
vtkm::cont::ArrayHandle<vtkm::Particle>& seeds,
const RK4Type& rk4,
vtkm::Id maxSteps,
vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>& result) const
{
vtkm::worklet::ParticleAdvection Worklet;
result = Worklet.Run(rk4, seeds, maxSteps);
}
//-----
// Specialization for Streamline worklet
template <>
template <>
inline void DataSetIntegratorBase<GridEvalType>::DoAdvect(
vtkm::cont::ArrayHandle<vtkm::Particle>& seeds,
const RK4Type& rk4,
vtkm::Id maxSteps,
vtkm::worklet::StreamlineResult<vtkm::Particle>& result) const
{
vtkm::worklet::Streamline Worklet;
result = Worklet.Run(rk4, seeds, maxSteps);
}
//-----
// Specialization for Pathline worklet
template <>
template <>
inline void DataSetIntegratorBase<TemporalGridEvalType>::DoAdvect(
vtkm::cont::ArrayHandle<vtkm::Particle>& seeds,
const TemporalRK4Type& rk4,
vtkm::Id maxSteps,
vtkm::worklet::StreamlineResult<vtkm::Particle>& result) const
{
vtkm::worklet::Streamline Worklet;
result = Worklet.Run(rk4, seeds, maxSteps);
}
}
}
} // namespace vtkm::filter::particleadvection
#endif

23
vtkm/filter/particleadvection/ParticleAdvectionAlgorithm.h
View File

@ -22,28 +22,29 @@ namespace filter
namespace particleadvection
{
class VTKM_ALWAYS_EXPORT ParticleAdvectionAlgorithm
: public AdvectorBaseAlgorithm<vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>>
{
using DataSetIntegratorType = vtkm::filter::particleadvection::DataSetIntegrator;
using DSIType = vtkm::filter::particleadvection::DataSetIntegrator;
class VTKM_ALWAYS_EXPORT ParticleAdvectionAlgorithm
: public AdvectorBaseAlgorithm<DSIType, vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>>
{
public:
ParticleAdvectionAlgorithm(const vtkm::filter::particleadvection::BoundsMap& bm,
const std::vector<DataSetIntegratorType>& blocks)
: AdvectorBaseAlgorithm<vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>>(bm, blocks)
const std::vector<DSIType>& blocks)
: AdvectorBaseAlgorithm<DSIType, vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>>(bm,
blocks)
{
}
};
class VTKM_ALWAYS_EXPORT ParticleAdvectionThreadedAlgorithm
: public AdvectorBaseThreadedAlgorithm<vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>>
: public AdvectorBaseThreadedAlgorithm<DSIType,
vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>>
{
using DataSetIntegratorType = vtkm::filter::particleadvection::DataSetIntegrator;
public:
ParticleAdvectionThreadedAlgorithm(const vtkm::filter::particleadvection::BoundsMap& bm,
const std::vector<DataSetIntegratorType>& blocks)
: AdvectorBaseThreadedAlgorithm<vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>>(bm,
const std::vector<DSIType>& blocks)
: AdvectorBaseThreadedAlgorithm<DSIType,
vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>>(bm,
blocks)
{
}

49
vtkm/filter/particleadvection/StreamlineAlgorithm.h
View File

@ -22,28 +22,55 @@ namespace filter
namespace particleadvection
{
class VTKM_ALWAYS_EXPORT StreamlineAlgorithm
: public AdvectorBaseAlgorithm<vtkm::worklet::StreamlineResult<vtkm::Particle>>
{
using DataSetIntegratorType = vtkm::filter::particleadvection::DataSetIntegrator;
using DSIType = vtkm::filter::particleadvection::DataSetIntegrator;
using TDSIType = vtkm::filter::particleadvection::TemporalDataSetIntegrator;
class VTKM_ALWAYS_EXPORT StreamlineAlgorithm
: public AdvectorBaseAlgorithm<DSIType, vtkm::worklet::StreamlineResult<vtkm::Particle>>
{
public:
StreamlineAlgorithm(const vtkm::filter::particleadvection::BoundsMap& bm,
const std::vector<DataSetIntegratorType>& blocks)
: AdvectorBaseAlgorithm<vtkm::worklet::StreamlineResult<vtkm::Particle>>(bm, blocks)
const std::vector<DSIType>& blocks)
: AdvectorBaseAlgorithm<DSIType, vtkm::worklet::StreamlineResult<vtkm::Particle>>(bm, blocks)
{
}
};
class VTKM_ALWAYS_EXPORT StreamlineThreadedAlgorithm
: public AdvectorBaseThreadedAlgorithm<vtkm::worklet::StreamlineResult<vtkm::Particle>>
: public AdvectorBaseThreadedAlgorithm<DSIType, vtkm::worklet::StreamlineResult<vtkm::Particle>>
{
using DataSetIntegratorType = vtkm::filter::particleadvection::DataSetIntegrator;
public:
StreamlineThreadedAlgorithm(const vtkm::filter::particleadvection::BoundsMap& bm,
const std::vector<DataSetIntegratorType>& blocks)
: AdvectorBaseThreadedAlgorithm<vtkm::worklet::StreamlineResult<vtkm::Particle>>(bm, blocks)
const std::vector<DSIType>& blocks)
: AdvectorBaseThreadedAlgorithm<DSIType, vtkm::worklet::StreamlineResult<vtkm::Particle>>(
bm,
blocks)
{
}
};
//pathline
class VTKM_ALWAYS_EXPORT PathlineAlgorithm
: public AdvectorBaseAlgorithm<TDSIType, vtkm::worklet::StreamlineResult<vtkm::Particle>>
{
public:
PathlineAlgorithm(const vtkm::filter::particleadvection::BoundsMap& bm,
const std::vector<TDSIType>& blocks)
: AdvectorBaseAlgorithm<TDSIType, vtkm::worklet::StreamlineResult<vtkm::Particle>>(bm, blocks)
{
}
};
class VTKM_ALWAYS_EXPORT PathlineThreadedAlgorithm
: public AdvectorBaseThreadedAlgorithm<TDSIType, vtkm::worklet::StreamlineResult<vtkm::Particle>>
{
public:
PathlineThreadedAlgorithm(const vtkm::filter::particleadvection::BoundsMap& bm,
const std::vector<TDSIType>& blocks)
: AdvectorBaseThreadedAlgorithm<TDSIType, vtkm::worklet::StreamlineResult<vtkm::Particle>>(
bm,
blocks)
{
}
};

1
vtkm/filter/testing/UnitTestStreamlineFilter.cxx
View File

@ -554,6 +554,7 @@ void TestStreamlineFilters()
TestStreamline();
TestPathlineSimple();
TestPathline();
for (auto useSL : flags)
TestAMRStreamline(useSL);

222
vtkm/filter/testing/UnitTestStreamlineFilterMPI.cxx
View File

@ -10,14 +10,21 @@
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/ParticleAdvection.h>
#include <vtkm/filter/Pathline.h>
#include <vtkm/filter/Streamline.h>
#include <vtkm/worklet/testing/GenerateTestDataSets.h>
#include <vtkm/thirdparty/diy/diy.h>
#include <vtkm/worklet/testing/GenerateTestDataSets.h>
namespace
{
enum FilterType
{
PARTICLE_ADVECTION,
STREAMLINE,
PATHLINE
};
vtkm::cont::ArrayHandle<vtkm::Vec3f> CreateConstantVectorField(vtkm::Id num, const vtkm::Vec3f& vec)
{
vtkm::cont::ArrayHandleConstant<vtkm::Vec3f> vecConst;
@ -36,7 +43,22 @@ void AddVectorFields(vtkm::cont::PartitionedDataSet& pds,
ds.AddPointField(fieldName, CreateConstantVectorField(ds.GetNumberOfPoints(), vec));
}
void TestAMRStreamline(bool useSL, bool useThreaded)
template <typename FilterType>
void SetFilter(FilterType& filter,
vtkm::FloatDefault stepSize,
vtkm::Id numSteps,
const std::string& fieldName,
vtkm::cont::ArrayHandle<vtkm::Particle> seedArray,
bool useThreaded)
{
filter.SetStepSize(stepSize);
filter.SetNumberOfSteps(numSteps);
filter.SetSeeds(seedArray);
filter.SetActiveField(fieldName);
filter.SetUseThreadedAlgorithm(useThreaded);
}
void TestAMRStreamline(FilterType fType, bool useThreaded)
{
auto comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
if (comm.size() < 2)
@ -102,16 +124,30 @@ void TestAMRStreamline(bool useSL, bool useThreaded)
seedArray = vtkm::cont::make_ArrayHandle({ p0, p1 });
vtkm::Id numSeeds = seedArray.GetNumberOfValues();
if (useSL)
{
vtkm::filter::Streamline filter;
filter.SetUseThreadedAlgorithm(useThreaded);
filter.SetStepSize(0.1f);
filter.SetNumberOfSteps(100000);
filter.SetSeeds(seedArray);
vtkm::FloatDefault stepSize = 0.1f;
vtkm::Id numSteps = 100000;
filter.SetActiveField(fieldName);
auto out = filter.Execute(pds);
if (fType == STREAMLINE || fType == PATHLINE)
{
vtkm::cont::PartitionedDataSet out;
if (fType == STREAMLINE)
{
vtkm::filter::Streamline streamline;
SetFilter(streamline, stepSize, numSteps, fieldName, seedArray, useThreaded);
out = streamline.Execute(pds);
}
else if (fType == PATHLINE)
{
vtkm::filter::Pathline pathline;
SetFilter(pathline, stepSize, numSteps, fieldName, seedArray, useThreaded);
//Create timestep 2
auto pds2 = vtkm::cont::PartitionedDataSet(pds);
pathline.SetPreviousTime(0);
pathline.SetNextTime(10);
pathline.SetNextDataSet(pds2);
out = pathline.Execute(pds);
}
if (comm.rank() <= 1)
VTKM_TEST_ASSERT(out.GetNumberOfPartitions() == 1, "Wrong number of partitions in output");
@ -185,7 +221,7 @@ void TestAMRStreamline(bool useSL, bool useThreaded)
"Seed final location in wrong location");
}
}
else
else if (fType == PARTICLE_ADVECTION)
{
vtkm::filter::ParticleAdvection filter;
filter.SetUseThreadedAlgorithm(useThreaded);
@ -221,7 +257,44 @@ void TestAMRStreamline(bool useSL, bool useThreaded)
}
}
void TestPartitionedDataSet(vtkm::Id nPerRank, bool useGhost, bool useSL, bool useThreaded)
void ValidateOutput(const vtkm::cont::DataSet& out,
vtkm::Id numSeeds,
vtkm::Range xMaxRange,
FilterType fType)
{
//Validate the result is correct.
VTKM_TEST_ASSERT(out.GetNumberOfCoordinateSystems() == 1,
"Wrong number of coordinate systems in the output dataset");
vtkm::cont::DynamicCellSet dcells = out.GetCellSet();
VTKM_TEST_ASSERT(dcells.GetNumberOfCells() == numSeeds, "Wrong number of cells");
auto coords = out.GetCoordinateSystem().GetDataAsMultiplexer();
auto ptPortal = coords.ReadPortal();
if (fType == STREAMLINE || fType == PATHLINE)
{
vtkm::cont::CellSetExplicit<> explicitCells;
VTKM_TEST_ASSERT(dcells.IsType<vtkm::cont::CellSetExplicit<>>(), "Wrong cell type.");
explicitCells = dcells.Cast<vtkm::cont::CellSetExplicit<>>();
for (vtkm::Id j = 0; j < numSeeds; j++)
{
vtkm::cont::ArrayHandle<vtkm::Id> indices;
explicitCells.GetIndices(j, indices);
vtkm::Id nPts = indices.GetNumberOfValues();
auto iPortal = indices.ReadPortal();
vtkm::Vec3f lastPt = ptPortal.Get(iPortal.Get(nPts - 1));
VTKM_TEST_ASSERT(xMaxRange.Contains(lastPt[0]), "Wrong end point for seed");
}
}
else if (fType == PARTICLE_ADVECTION)
{
VTKM_TEST_ASSERT(out.GetNumberOfPoints() == numSeeds, "Wrong number of coordinates");
for (vtkm::Id i = 0; i < numSeeds; i++)
VTKM_TEST_ASSERT(xMaxRange.Contains(ptPortal.Get(i)[0]), "Wrong end point for seed");
}
}
void TestPartitionedDataSet(vtkm::Id nPerRank, bool useGhost, FilterType fType, bool useThreaded)
{
auto comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
@ -243,6 +316,9 @@ void TestPartitionedDataSet(vtkm::Id nPerRank, bool useGhost, bool useSL, bool u
z1 = z1 + 1;
}
vtkm::FloatDefault time0 = 0;
vtkm::FloatDefault time1 = (numDims - 1) * (nPerRank * comm.size());
std::vector<vtkm::Bounds> bounds;
for (vtkm::Id i = 0; i < nPerRank; i++)
{
@ -251,14 +327,28 @@ void TestPartitionedDataSet(vtkm::Id nPerRank, bool useGhost, bool useSL, bool u
x1 += dx;
}
std::vector<vtkm::cont::PartitionedDataSet> allPDs;
std::vector<vtkm::Range> xMaxRanges;
for (const auto& b : bounds)
{
auto xMax = b.X.Max;
if (useGhost)
xMax = xMax - 1;
xMaxRanges.push_back(vtkm::Range(xMax, xMax + static_cast<vtkm::FloatDefault>(.5)));
}
std::vector<vtkm::cont::PartitionedDataSet> allPDS, allPDS2;
const vtkm::Id3 dims(numDims, numDims, numDims);
allPDs = vtkm::worklet::testing::CreateAllDataSets(bounds, dims, useGhost);
allPDS = vtkm::worklet::testing::CreateAllDataSets(bounds, dims, useGhost);
allPDS2 = vtkm::worklet::testing::CreateAllDataSets(bounds, dims, useGhost);
vtkm::Vec3f vecX(1, 0, 0);
std::string fieldName = "vec";
for (auto& pds : allPDs)
vtkm::FloatDefault stepSize = 0.1f;
vtkm::Id numSteps = 100000;
for (std::size_t n = 0; n < allPDS.size(); n++)
{
auto pds = allPDS[n];
AddVectorFields(pds, fieldName, vecX);
vtkm::cont::ArrayHandle<vtkm::Particle> seedArray;
@ -266,107 +356,65 @@ void TestPartitionedDataSet(vtkm::Id nPerRank, bool useGhost, bool useSL, bool u
vtkm::Particle(vtkm::Vec3f(.2f, 2.0f, .2f), 1) });
vtkm::Id numSeeds = seedArray.GetNumberOfValues();
if (useSL)
if (fType == STREAMLINE)
{
vtkm::filter::Streamline streamline;
streamline.SetUseThreadedAlgorithm(useThreaded);
streamline.SetStepSize(0.1f);
streamline.SetNumberOfSteps(100000);
streamline.SetSeeds(seedArray);
streamline.SetActiveField(fieldName);
SetFilter(streamline, stepSize, numSteps, fieldName, seedArray, useThreaded);
auto out = streamline.Execute(pds);
for (vtkm::Id i = 0; i < nPerRank; i++)
{
auto inputDS = pds.GetPartition(i);
auto outputDS = out.GetPartition(i);
VTKM_TEST_ASSERT(outputDS.GetNumberOfCoordinateSystems() == 1,
"Wrong number of coordinate systems in the output dataset");
vtkm::cont::DynamicCellSet dcells = outputDS.GetCellSet();
VTKM_TEST_ASSERT(dcells.GetNumberOfCells() == numSeeds, "Wrong number of cells");
auto coords = outputDS.GetCoordinateSystem().GetDataAsMultiplexer();
auto ptPortal = coords.ReadPortal();
vtkm::cont::CellSetExplicit<> explicitCells;
VTKM_TEST_ASSERT(dcells.IsType<vtkm::cont::CellSetExplicit<>>(), "Wrong cell type.");
explicitCells = dcells.Cast<vtkm::cont::CellSetExplicit<>>();
vtkm::FloatDefault xMax = bounds[static_cast<std::size_t>(i)].X.Max;
if (useGhost)
xMax = xMax - 1;
vtkm::Range xMaxRange(xMax, xMax + static_cast<vtkm::FloatDefault>(.5));
for (vtkm::Id j = 0; j < numSeeds; j++)
{
vtkm::cont::ArrayHandle<vtkm::Id> indices;
explicitCells.GetIndices(j, indices);
vtkm::Id nPts = indices.GetNumberOfValues();
auto iPortal = indices.ReadPortal();
vtkm::Vec3f lastPt = ptPortal.Get(iPortal.Get(nPts - 1));
VTKM_TEST_ASSERT(xMaxRange.Contains(lastPt[0]), "Wrong end point for seed");
}
}
ValidateOutput(out.GetPartition(i), numSeeds, xMaxRanges[i], fType);
}
else
else if (fType == PARTICLE_ADVECTION)
{
vtkm::filter::ParticleAdvection particleAdvection;
particleAdvection.SetUseThreadedAlgorithm(useThreaded);
particleAdvection.SetStepSize(0.1f);
particleAdvection.SetNumberOfSteps(100000);
particleAdvection.SetSeeds(seedArray);
particleAdvection.SetActiveField(fieldName);
SetFilter(particleAdvection, stepSize, numSteps, fieldName, seedArray, useThreaded);
auto out = particleAdvection.Execute(pds);
//Particles end up in last rank.
if (comm.rank() == comm.size() - 1)
{
VTKM_TEST_ASSERT(out.GetNumberOfPartitions() == 1, "Wrong number of partitions in output");
auto ds = out.GetPartition(0);
//Validate the result is correct.
VTKM_TEST_ASSERT(ds.GetNumberOfCoordinateSystems() == 1,
"Wrong number of coordinate systems in the output dataset");
vtkm::FloatDefault xMax = bounds[bounds.size() - 1].X.Max;
if (useGhost)
xMax = xMax - 1;
vtkm::Range xMaxRange(xMax, xMax + static_cast<vtkm::FloatDefault>(.5));
auto coords = ds.GetCoordinateSystem().GetDataAsMultiplexer();
VTKM_TEST_ASSERT(ds.GetNumberOfPoints() == numSeeds, "Wrong number of coordinates");
auto ptPortal = coords.ReadPortal();
for (vtkm::Id i = 0; i < numSeeds; i++)
VTKM_TEST_ASSERT(xMaxRange.Contains(ptPortal.Get(i)[0]), "Wrong end point for seed");
vtkm::cont::DynamicCellSet dcells = ds.GetCellSet();
VTKM_TEST_ASSERT(dcells.GetNumberOfCells() == numSeeds, "Wrong number of cells");
ValidateOutput(out.GetPartition(0), numSeeds, xMaxRanges[xMaxRanges.size() - 1], fType);
}
else
VTKM_TEST_ASSERT(out.GetNumberOfPartitions() == 0, "Wrong number of partitions in output");
}
else if (fType == PATHLINE)
{
auto pds2 = allPDS2[n];
AddVectorFields(pds2, fieldName, vecX);
vtkm::filter::Pathline pathline;
SetFilter(pathline, stepSize, numSteps, fieldName, seedArray, useThreaded);
pathline.SetPreviousTime(time0);
pathline.SetNextTime(time1);
pathline.SetNextDataSet(pds2);
auto out = pathline.Execute(pds);
for (vtkm::Id i = 0; i < nPerRank; i++)
ValidateOutput(out.GetPartition(i), numSeeds, xMaxRanges[i], fType);
}
}
}
void TestStreamlineFiltersMPI()
{
std::vector<bool> flags = { true, false };
std::vector<FilterType> filterTypes = { PARTICLE_ADVECTION, STREAMLINE, PATHLINE };
for (int n = 1; n < 3; n++)
{
for (auto useGhost : flags)
for (auto useSL : flags)
for (auto fType : filterTypes)
for (auto useThreaded : flags)
TestPartitionedDataSet(n, useGhost, useSL, useThreaded);
TestPartitionedDataSet(n, useGhost, fType, useThreaded);
}
for (auto useSL : flags)
for (auto fType : filterTypes)
for (auto useThreaded : flags)
TestAMRStreamline(useSL, useThreaded);
TestAMRStreamline(fType, useThreaded);
}
}

15
vtkm/worklet/particleadvection/TemporalGridEvaluators.h
View File

@ -121,6 +121,21 @@ private:
public:
VTKM_CONT TemporalGridEvaluator() = default;
VTKM_CONT TemporalGridEvaluator(const vtkm::cont::DataSet& ds1,
const vtkm::FloatDefault t1,
const FieldType& field1,
const vtkm::cont::DataSet& ds2,
const vtkm::FloatDefault t2,
const FieldType& field2)
: EvaluatorOne(GridEvaluator(ds1, field1))
, EvaluatorTwo(GridEvaluator(ds2, field2))
, TimeOne(t1)
, TimeTwo(t2)
{
}
VTKM_CONT TemporalGridEvaluator(GridEvaluator& evaluatorOne,
const vtkm::FloatDefault timeOne,
GridEvaluator& evaluatorTwo,