Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-16 17:22:55 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Added a version that directly uses ArrayHandle<vtkm::Particle>

Browse Source
This commit is contained in:
Dave Pugmire 2019-09-18 14:08:26 -04:00
parent 0c70950596
commit 28a73494b1
5 changed files with 183 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
examples/particle_advection/ParticleAdvection.cxx
View File

@ -82,7 +82,7 @@ int main(int argc, char** argv)
//create seeds randomly placed withing the bounding box of the data.
vtkm::Bounds bounds = ds.GetCoordinateSystem().GetBounds();
std::vector<vtkm::Vec3f> seeds;
std::vector<vtkm::Particle> particles;
std::vector<vtkm::Particle> particles, particles2;
for (int i = 0; i < numSeeds; i++)
{
@ -101,10 +101,12 @@ int main(int argc, char** argv)
pa.Pos = p;
pa.NumSteps = 0;
particles.push_back(pa);
particles2.push_back(pa);
}
auto seedArrayPts = vtkm::cont::make_ArrayHandle(seeds);
auto seedArrayPar = vtkm::cont::make_ArrayHandle(particles);
auto seedArrayPar2 = vtkm::cont::make_ArrayHandle(particles2);
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldHandle>;
@ -131,12 +133,21 @@ int main(int argc, char** argv)
if (1)
{
auto s = std::chrono::system_clock::now();
auto res1 = paPART.Run(rk4, seedArrayPar, numSteps);
auto res1 = paPART.Run(rk4, seedArrayPar, numSteps, 0);
auto e = std::chrono::system_clock::now();
dT = e - s;
std::cout << "AOS_Time= " << dT.count() << std::endl;
}
if (1)
{
auto s = std::chrono::system_clock::now();
auto res1 = paPART.Run(rk4, seedArrayPar2, numSteps, 1);
auto e = std::chrono::system_clock::now();
dT = e - s;
std::cout << "AOS2_Time= " << dT.count() << std::endl;
}
/*
//compute streamlines
vtkm::filter::Streamline streamline;

31
vtkm/Particle.h
View File

@ -29,12 +29,37 @@ enum ParticleStatus : vtkm::Id
class Particle
{
public:
Particle() {}
Particle(const vtkm::Vec3f& p, vtkm::Id id, vtkm::Id numSteps)
VTKM_EXEC_CONT
Particle()
: Pos()
, ID(-1)
, NumSteps(0)
, Status(vtkm::ParticleStatus::UNDEFINED)
, Time(0)
{
}
VTKM_EXEC_CONT
Particle(const vtkm::Particle& p)
: Pos(p.Pos)
, ID(p.ID)
, NumSteps(p.NumSteps)
, Status(p.Status)
, Time(p.Time)
{
}
VTKM_EXEC_CONT
Particle(const vtkm::Vec3f& p,
vtkm::Id id,
vtkm::Id numSteps,
vtkm::Id status = vtkm::ParticleStatus::UNDEFINED,
vtkm::FloatDefault time = 0)
: Pos(p)
, ID(id)
, NumSteps(numSteps)
, Status(vtkm::ParticleStatus::UNDEFINED)
, Status(status)
, Time(time)
{
}

24
vtkm/worklet/ParticleAdvection.h
View File

@ -139,27 +139,17 @@ public:
template <typename IntegratorType, typename ParticleStorage>
ParticleAdvectionResult Run(const IntegratorType& it,
vtkm::cont::ArrayHandle<vtkm::Particle, ParticleStorage>& particles,
const vtkm::Id& MaxSteps)
const vtkm::Id& MaxSteps,
int which)
{
vtkm::worklet::particleadvection::ParticleAdvectionWorklet<IntegratorType> worklet;
worklet.Run(it, particles, MaxSteps);
if (which == 0)
worklet.Run(it, particles, MaxSteps);
else
worklet.Run2(it, particles, MaxSteps);
return ParticleAdvectionResult();
/*
vtkm::Id numSeeds = static_cast<vtkm::Id>(particles.GetNumberOfValues());
vtkm::cont::ArrayHandle<vtkm::Id> status;
//Allocate status arrays.
vtkm::cont::ArrayHandleConstant<vtkm::Id> statusOK(static_cast<vtkm::Id>(1), numSeeds);
status.Allocate(numSeeds);
vtkm::cont::ArrayCopy(statusOK, status);
worklet.Run(it, pts, nSteps, status, inputSteps, inputTime);
//Create output.
return ParticleAdvectionResult(pts, status, inputSteps, inputTime);
*/
}
};

111
vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h
View File

@ -101,6 +101,92 @@ public:
};
class ParticleAdvectWorkletAOS2 : public vtkm::worklet::WorkletMapField
{
public:
using ControlSignature = void(ExecObject integrator, FieldInOut particles, FieldIn maxSteps);
using ExecutionSignature = void(_1, _2, _3);
using InputDomain = _3;
template <typename IntegratorType, typename ParticleType>
VTKM_EXEC void operator()(const IntegratorType* integrator,
ParticleType& particle,
const vtkm::Id& maxSteps) const
{
vtkm::Vec3f inpos = particle.Pos, outpos;
vtkm::FloatDefault time = particle.Time;
IntegratorStatus status;
bool tookAnySteps = false;
while ((particle.Status & ParticleStatus::SUCCESS) != 0)
{
status = integrator->Step(inpos, time, outpos);
// If the status is OK, we only need to check if the particle
// has completed the maximum steps required.
if (status == IntegratorStatus::SUCCESS)
{
//particle.TakeStep(idx, outpos, time);
particle.Pos = outpos;
particle.NumSteps++;
particle.Time = time;
if (particle.NumSteps == maxSteps)
{
particle.Status &= ~ParticleStatus::SUCCESS;
particle.Status |= ParticleStatus::TERMINATED;
}
// This is to keep track of the particle's time.
// This is what the Evaluator uses to determine if the particle
// has exited temporal boundary.
inpos = outpos;
tookAnySteps = true;
}
// If the particle is at spatial or temporal boundary, take steps to just
// push it a little out of the boundary so that it will start advection in
// another domain, or in another time slice. Taking small steps enables
// reducing the error introduced at spatial or temporal boundaries.
else if (status == IntegratorStatus::OUTSIDE_TEMPORAL_BOUNDS)
{
particle.Status &= ~ParticleStatus::SUCCESS;
particle.Status |= ParticleStatus::EXIT_TEMPORAL_BOUNDARY;
break;
}
else if (status == IntegratorStatus::OUTSIDE_SPATIAL_BOUNDS)
{
status = integrator->SmallStep(inpos, time, outpos);
//particle.TakeStep(idx, outpos, time);
particle.Pos = outpos;
particle.NumSteps++;
particle.Time = time;
if (status == IntegratorStatus::OUTSIDE_TEMPORAL_BOUNDS)
{
particle.Status &= ~ParticleStatus::SUCCESS;
particle.Status |= ParticleStatus::EXIT_TEMPORAL_BOUNDARY;
break;
}
else if (status == IntegratorStatus::OUTSIDE_SPATIAL_BOUNDS)
{
particle.Status &= ~ParticleStatus::SUCCESS;
particle.Status |= ParticleStatus::EXIT_SPATIAL_BOUNDARY;
break;
}
else if (status == IntegratorStatus::FAIL)
{
particle.Status &= ~ParticleStatus::SUCCESS;
particle.Status |= ParticleStatus::FAIL;
break;
}
}
}
if (tookAnySteps)
particle.Status |= ParticleStatus::TOOK_ANY_STEPS;
else
particle.Status &= ~ParticleStatus::TOOK_ANY_STEPS;
}
};
class ParticleAdvectWorklet : public vtkm::worklet::WorkletMapField
{
public:
@ -228,9 +314,32 @@ public:
//Invoke particle advection worklet
ParticleWorkletDispatchType particleWorkletDispatch;
//problem
particleWorkletDispatch.Invoke(idxArray, integrator, particles);
}
//AOS version
void Run2(const IntegratorType& integrator,
vtkm::cont::ArrayHandle<vtkm::Particle>& particles,
const vtkm::Id& MaxSteps)
{
using ParticleAdvectWorkletType = vtkm::worklet::particleadvection::ParticleAdvectWorkletAOS2;
using ParticleWorkletDispatchType =
typename vtkm::worklet::DispatcherMapField<ParticleAdvectWorkletType>;
vtkm::Id numSeeds = static_cast<vtkm::Id>(particles.GetNumberOfValues());
//Create and invoke the particle advection.
vtkm::cont::ArrayHandleConstant<vtkm::Id> maxSteps(MaxSteps, numSeeds);
#ifdef VTKM_CUDA
// This worklet needs some extra space on CUDA.
vtkm::cont::cuda::ScopedCudaStackSize stack(16 * 1024);
(void)stack;
#endif // VTKM_CUDA
//Invoke particle advection worklet
ParticleWorkletDispatchType particleWorkletDispatch;
particleWorkletDispatch.Invoke(integrator, particles, maxSteps);
}
};
namespace detail

49
vtkm/worklet/testing/UnitTestParticleAdvection.cxx
View File

@ -658,7 +658,7 @@ void TestParticleStatus()
VTKM_TEST_ASSERT(tookStep1 == 0, "Particle took a step when it should not have.");
}
#if 0
double TestParticleAOS(bool doAOS)
{
vtkm::Bounds bounds(0, 1, 0, 1, 0, 1);
@ -696,37 +696,38 @@ double TestParticleAOS(bool doAOS)
int n = 10000;
for (int i = 0; i < n; i++)
{
vtkm::Particle p;
p.ID = i;
p.Status = vtkm::ParticleStatus::SUCCESS;
auto pt = RandomPoint(bounds);
p.Pos = pt;
particles.push_back(p);
pts.push_back(pt);
vtkm::Particle p;
p.ID = i;
p.Status = vtkm::ParticleStatus::SUCCESS;
auto pt = RandomPoint(bounds);
p.Pos = pt;
particles.push_back(p);
pts.push_back(pt);
}
std::chrono::duration<double> dT;
if (doAOS)
{
auto seedsArray = vtkm::cont::make_ArrayHandle(particles);
//vtkm::cont::printSummary_ArrayHandle(seedsArray, std::cout, true);
auto s = std::chrono::system_clock::now();
auto res = pa.Run(rk4, seedsArray, maxSteps);
auto e = std::chrono::system_clock::now();
dT = e - s;
auto seedsArray = vtkm::cont::make_ArrayHandle(particles);
//vtkm::cont::printSummary_ArrayHandle(seedsArray, std::cout, true);
auto s = std::chrono::system_clock::now();
auto res = pa.Run(rk4, seedsArray, maxSteps);
auto e = std::chrono::system_clock::now();
dT = e-s;
}
else
{
auto seedsArray = vtkm::cont::make_ArrayHandle(pts);
//vtkm::cont::printSummary_ArrayHandle(seedsArray, std::cout, true);
auto s = std::chrono::system_clock::now();
auto res = pa.Run(rk4, seedsArray, maxSteps);
auto e = std::chrono::system_clock::now();
dT = e - s;
auto seedsArray = vtkm::cont::make_ArrayHandle(pts);
//vtkm::cont::printSummary_ArrayHandle(seedsArray, std::cout, true);
auto s = std::chrono::system_clock::now();
auto res = pa.Run(rk4, seedsArray, maxSteps);
auto e = std::chrono::system_clock::now();
dT = e-s;
}
return dT.count();
}
#endif
void TestParticleAdvection()
{
@ -734,11 +735,11 @@ void TestParticleAdvection()
// TestParticleWorklets();
// TestParticleStatus();
double aosTime = TestParticleAOS(true);
double soaTime = TestParticleAOS(false);
// double aosTime = TestParticleAOS(true);
// double soaTime = TestParticleAOS(false);
std::cout << "AOS_time= " << aosTime << std::endl;
std::cout << "SOA_time= " << soaTime << std::endl;
// std::cout<<"AOS_time= "<< aosTime<<std::endl;
// std::cout<<"SOA_time= "<< soaTime<<std::endl;
}
int UnitTestParticleAdvection(int argc, char* argv[])