diff --git a/examples/particle_advection/ParticleAdvection.cxx b/examples/particle_advection/ParticleAdvection.cxx index af64d2e9d..bcd9d0475 100644 --- a/examples/particle_advection/ParticleAdvection.cxx +++ b/examples/particle_advection/ParticleAdvection.cxx @@ -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 seeds; - std::vector particles; + std::vector 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; using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator; @@ -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; diff --git a/vtkm/Particle.h b/vtkm/Particle.h index 7dcc76654..3c7bea46c 100644 --- a/vtkm/Particle.h +++ b/vtkm/Particle.h @@ -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) { } diff --git a/vtkm/worklet/ParticleAdvection.h b/vtkm/worklet/ParticleAdvection.h index 8f01aa7f8..d108580f9 100644 --- a/vtkm/worklet/ParticleAdvection.h +++ b/vtkm/worklet/ParticleAdvection.h @@ -139,27 +139,17 @@ public: template ParticleAdvectionResult Run(const IntegratorType& it, vtkm::cont::ArrayHandle& particles, - const vtkm::Id& MaxSteps) + const vtkm::Id& MaxSteps, + int which) { vtkm::worklet::particleadvection::ParticleAdvectionWorklet 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(particles.GetNumberOfValues()); - - vtkm::cont::ArrayHandle status; - //Allocate status arrays. - vtkm::cont::ArrayHandleConstant statusOK(static_cast(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); -*/ } }; diff --git a/vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h b/vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h index 510b2f608..d4d727bfe 100644 --- a/vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h +++ b/vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h @@ -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 + 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& particles, + const vtkm::Id& MaxSteps) + { + using ParticleAdvectWorkletType = vtkm::worklet::particleadvection::ParticleAdvectWorkletAOS2; + using ParticleWorkletDispatchType = + typename vtkm::worklet::DispatcherMapField; + + vtkm::Id numSeeds = static_cast(particles.GetNumberOfValues()); + //Create and invoke the particle advection. + vtkm::cont::ArrayHandleConstant 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 diff --git a/vtkm/worklet/testing/UnitTestParticleAdvection.cxx b/vtkm/worklet/testing/UnitTestParticleAdvection.cxx index 1fa9514c6..fa5e084bb 100644 --- a/vtkm/worklet/testing/UnitTestParticleAdvection.cxx +++ b/vtkm/worklet/testing/UnitTestParticleAdvection.cxx @@ -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 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<