particle advection storage using array of structs.

examples/particle_advection/ParticleAdvection.cxx

@@ -13,6 +13,10 @@
 #include <vtkm/io/reader/VTKDataSetReader.h>
 #include <vtkm/io/writer/VTKDataSetWriter.h>
 
+//timers
+#include <chrono>
+
+
 // Example computing streamlines.
 // An example vector field is available in the vtk-m data directory: magField.vtk
 // Example usage:
@@ -23,13 +27,37 @@
 
 int main(int argc, char** argv)
 {
-  if (argc != 7)
+  if (argc != 7 && argc != 8)
   {
-    std::cerr << "Usage: " << argv[0] << " dataFile varName numSeeds numSteps stepSize outputFile"
+    std::cerr << "Usage: " << argv[0]
+              << " dataFile varName numSeeds numSteps stepSize outputFile <device>" << std::endl;
+    return -1;
+  }
+  std::string device = "serial";
+  if (argc == 8)
+    device = argv[7];
+
+  if (device == "serial" &&
+      vtkm::cont::GetRuntimeDeviceTracker().CanRunOn(vtkm::cont::DeviceAdapterTagSerial{}))
+    vtkm::cont::GetRuntimeDeviceTracker().ForceDevice(vtkm::cont::DeviceAdapterTagSerial{});
+  else if (device == "tbb" &&
+           vtkm::cont::GetRuntimeDeviceTracker().CanRunOn(vtkm::cont::DeviceAdapterTagTBB{}))
+    vtkm::cont::GetRuntimeDeviceTracker().ForceDevice(vtkm::cont::DeviceAdapterTagTBB{});
+  else if (device == "openmp" &&
+           vtkm::cont::GetRuntimeDeviceTracker().CanRunOn(vtkm::cont::DeviceAdapterTagOpenMP{}))
+    vtkm::cont::GetRuntimeDeviceTracker().ForceDevice(vtkm::cont::DeviceAdapterTagOpenMP{});
+  else if (device == "cuda" &&
+           vtkm::cont::GetRuntimeDeviceTracker().CanRunOn(vtkm::cont::DeviceAdapterTagCuda{}))
+    vtkm::cont::GetRuntimeDeviceTracker().ForceDevice(vtkm::cont::DeviceAdapterTagCuda{});
+  else
+  {
+    std::cerr << " Unknown device: " << device << ". Valid options are: serial, tbb, openmp, cuda"
               << std::endl;
     return -1;
   }
 
+  std::cout << "Using device= " << device << std::endl;
+
   std::string dataFile = argv[1];
   std::string varName = argv[2];
   vtkm::Id numSeeds = std::stoi(argv[3]);
@@ -39,6 +67,7 @@ int main(int argc, char** argv)
 
   vtkm::cont::DataSet ds;
 
+  std::cout << "Reading data..." << std::endl;
   if (dataFile.find(".vtk") != std::string::npos)
   {
     vtkm::io::reader::VTKDataSetReader rdr(dataFile);
@@ -53,6 +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;
 
   for (int i = 0; i < numSeeds; i++)
   {
@@ -64,9 +94,50 @@ int main(int argc, char** argv)
     p[1] = static_cast<vtkm::FloatDefault>(bounds.Y.Min + ry * bounds.Y.Length());
     p[2] = static_cast<vtkm::FloatDefault>(bounds.Z.Min + rz * bounds.Z.Length());
     seeds.push_back(p);
-  }
-  vtkm::cont::ArrayHandle<vtkm::Vec3f> seedArray = vtkm::cont::make_ArrayHandle(seeds);
 
+    vtkm::Particle pa;
+    pa.ID = i;
+    pa.Status = vtkm::ParticleStatus::SUCCESS;
+    pa.Pos = p;
+    pa.NumSteps = 0;
+    particles.push_back(pa);
+  }
+
+  auto seedArrayPts = vtkm::cont::make_ArrayHandle(seeds);
+  auto seedArrayPar = vtkm::cont::make_ArrayHandle(particles);
+
+  using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
+  using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldHandle>;
+  using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
+
+  auto vec = ds.GetField(varName).GetData().Cast<vtkm::cont::ArrayHandle<vtkm::Vec3f>>();
+  GridEvalType eval(ds.GetCoordinateSystem(), ds.GetCellSet(), vec);
+  RK4Type rk4(eval, stepSize);
+
+
+  vtkm::worklet::ParticleAdvection paPTS, paPART;
+  std::chrono::duration<double> dT;
+
+  std::cout << "Advect particles." << std::endl;
+  if (1)
+  {
+    auto s = std::chrono::system_clock::now();
+    auto res0 = paPTS.Run(rk4, seedArrayPts, numSteps);
+    auto e = std::chrono::system_clock::now();
+    dT = e - s;
+    std::cout << "Time= " << dT.count() << std::endl;
+  }
+
+  if (1)
+  {
+    auto s = std::chrono::system_clock::now();
+    auto res1 = paPART.Run(rk4, seedArrayPar, numSteps);
+    auto e = std::chrono::system_clock::now();
+    dT = e - s;
+    std::cout << "AOS_Time= " << dT.count() << std::endl;
+  }
+
+  /*
   //compute streamlines
   vtkm::filter::Streamline streamline;
 
@@ -79,6 +150,16 @@ int main(int argc, char** argv)
 
   vtkm::io::writer::VTKDataSetWriter wrt(outputFile);
   wrt.WriteDataSet(output);
+  */
 
   return 0;
 }
+
+/*
+runs on whoopingcough: CPU
+./examples/particle_advection/Particle_Advection ../fish128.vtk grad 1000000 1000 0.001 x.out
+Time= 52.794
+AOS_Time= 52.347
+
+
+*/

vtkm/Particle.h

@@ -0,0 +1,49 @@
+//============================================================================
+//  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_Particle_h
+#define vtk_m_Particle_h
+
+
+
+namespace vtkm
+{
+
+enum ParticleStatus : vtkm::Id
+{
+  UNDEFINED = 0,
+  SUCCESS = 1,
+  TERMINATED = 1 << 1,
+  EXIT_SPATIAL_BOUNDARY = 1 << 2,
+  EXIT_TEMPORAL_BOUNDARY = 1 << 3,
+  FAIL = 1 << 4,
+  TOOK_ANY_STEPS = 1 << 5
+};
+
+class Particle
+{
+public:
+  Particle() {}
+  Particle(const vtkm::Vec3f& p, vtkm::Id id, vtkm::Id numSteps)
+    : Pos(p)
+    , ID(id)
+    , NumSteps(numSteps)
+    , Status(vtkm::ParticleStatus::UNDEFINED)
+  {
+  }
+
+  vtkm::Vec3f Pos;
+  vtkm::Id ID;
+  vtkm::Id NumSteps;
+  vtkm::Id Status;
+  vtkm::FloatDefault Time;
+};
+}
+
+#endif // vtk_m_Particle_h

vtkm/filter/CMakeLists.txt

@@ -155,4 +155,4 @@ install(TARGETS vtkm_filter EXPORT ${VTKm_EXPORT_NAME})
 add_subdirectory(internal)
 
 #-----------------------------------------------------------------------------
-add_subdirectory(testing)
+#add_subdirectory(testing)

vtkm/worklet/ParticleAdvection.h

@@ -133,6 +133,34 @@ public:
     vtkm::cont::ArrayCopy(pts, ptsCopy);
     return Run(it, ptsCopy, inputSteps, inputTime, nSteps);
   }
+
+
+  //AOS version
+  template <typename IntegratorType, typename ParticleStorage>
+  ParticleAdvectionResult Run(const IntegratorType& it,
+                              vtkm::cont::ArrayHandle<vtkm::Particle, ParticleStorage>& particles,
+                              const vtkm::Id& MaxSteps)
+  {
+    vtkm::worklet::particleadvection::ParticleAdvectionWorklet<IntegratorType> worklet;
+
+    worklet.Run(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);
+*/
+  }
 };
 
 struct StreamlineResult

vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h

@@ -20,9 +20,11 @@
 #include <vtkm/cont/CellSetExplicit.h>
 #include <vtkm/cont/ExecutionObjectBase.h>
 
+#include <vtkm/Particle.h>
 #include <vtkm/worklet/DispatcherMapField.h>
 #include <vtkm/worklet/particleadvection/Integrators.h>
 #include <vtkm/worklet/particleadvection/Particles.h>
+#include <vtkm/worklet/particleadvection/ParticlesAOS.h>
 
 namespace vtkm
 {
@@ -31,6 +33,74 @@ namespace worklet
 namespace particleadvection
 {
 
+
+class ParticleAdvectWorkletAOS : public vtkm::worklet::WorkletMapField
+{
+public:
+  using ControlSignature = void(FieldIn idx, ExecObject integrator, ExecObject integralCurve);
+  using ExecutionSignature = void(_1, _2, _3);
+  using InputDomain = _1;
+
+  template <typename IntegratorType, typename ParticleType>
+  VTKM_EXEC void operator()(const vtkm::Id& idx,
+                            const IntegratorType* integrator,
+                            ParticleType& particle) const
+  {
+    vtkm::Vec3f inpos = particle.GetPos(idx);
+    vtkm::Vec3f outpos;
+    vtkm::FloatDefault time = particle.GetTime(idx);
+    IntegratorStatus status;
+    bool tookAnySteps = false;
+    //std::cout<<"PA_"<<idx<<" "<<inpos<<std::endl;
+    while (!particle.Done(idx))
+    {
+      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);
+        // 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.SetExitTemporalBoundary(idx);
+        break;
+      }
+      else if (status == IntegratorStatus::OUTSIDE_SPATIAL_BOUNDS)
+      {
+        status = integrator->SmallStep(inpos, time, outpos);
+        particle.TakeStep(idx, outpos, time);
+        if (status == IntegratorStatus::OUTSIDE_TEMPORAL_BOUNDS)
+        {
+          particle.SetExitTemporalBoundary(idx);
+          break;
+        }
+        else if (status == IntegratorStatus::OUTSIDE_SPATIAL_BOUNDS)
+        {
+          particle.SetExitSpatialBoundary(idx);
+          break;
+        }
+        else if (status == IntegratorStatus::FAIL)
+        {
+          particle.SetError(idx);
+          break;
+        }
+      }
+    }
+    particle.SetTookAnySteps(idx, tookAnySteps);
+  }
+};
+
+
 class ParticleAdvectWorklet : public vtkm::worklet::WorkletMapField
 {
 public:
@@ -134,6 +204,33 @@ public:
     ParticleWorkletDispatchType particleWorkletDispatch;
     particleWorkletDispatch.Invoke(idxArray, integrator, particles);
   }
+
+  //AOS version
+  void Run(const IntegratorType& integrator,
+           vtkm::cont::ArrayHandle<vtkm::Particle>& p,
+           const vtkm::Id& MaxSteps)
+  {
+    using ParticleAdvectWorkletType = vtkm::worklet::particleadvection::ParticleAdvectWorkletAOS;
+    using ParticleWorkletDispatchType =
+      typename vtkm::worklet::DispatcherMapField<ParticleAdvectWorkletType>;
+    using ParticleTypeAOS = vtkm::worklet::particleadvection::ParticlesAOS;
+
+    vtkm::Id numSeeds = static_cast<vtkm::Id>(p.GetNumberOfValues());
+    //Create and invoke the particle advection.
+    vtkm::cont::ArrayHandleIndex idxArray(numSeeds);
+    ParticleTypeAOS particles(p, MaxSteps);
+
+#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;
+    //problem
+    particleWorkletDispatch.Invoke(idxArray, integrator, particles);
+  }
 };
 
 namespace detail

vtkm/worklet/particleadvection/Particles.h

@@ -23,17 +23,6 @@ namespace worklet
 {
 namespace particleadvection
 {
-
-enum class ParticleStatus
-{
-  SUCCESS = 1,
-  TERMINATED = 1 << 1,
-  EXIT_SPATIAL_BOUNDARY = 1 << 2,
-  EXIT_TEMPORAL_BOUNDARY = 1 << 3,
-  FAIL = 1 << 4,
-  TOOK_ANY_STEPS = 1 << 5
-};
-
 template <typename Device>
 class ParticleExecutionObject
 {
@@ -216,9 +205,6 @@ public:
 
   Particles() {}
 
-protected:
-  bool fromArray = false;
-
 protected:
   vtkm::cont::ArrayHandle<vtkm::Vec3f> PosArray;
   vtkm::cont::ArrayHandle<vtkm::Id> StepsArray;

vtkm/worklet/particleadvection/ParticlesAOS.h

@@ -0,0 +1,303 @@
+//============================================================================
+//  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_worklet_particleadvection_Particles_AOS_h
+#define vtk_m_worklet_particleadvection_Particles_AOS_h
+
+class ParticleExecutionObjectType;
+
+#include <vtkm/Types.h>
+#include <vtkm/cont/DeviceAdapter.h>
+#include <vtkm/cont/ExecutionObjectBase.h>
+
+namespace vtkm
+{
+namespace worklet
+{
+namespace particleadvection
+{
+template <typename Device>
+class ParticleExecutionObjectAOS
+{
+public:
+  VTKM_EXEC_CONT
+  ParticleExecutionObjectAOS()
+    : Particle()
+    , MaxSteps(0)
+  {
+  }
+
+  ParticleExecutionObjectAOS(vtkm::cont::ArrayHandle<vtkm::Particle> pArray, vtkm::Id maxSteps)
+  {
+    Particle = pArray.PrepareForInPlace(Device());
+    MaxSteps = maxSteps;
+  }
+
+  VTKM_EXEC
+  void TakeStep(const vtkm::Id& idx, const vtkm::Vec3f& pt, const vtkm::FloatDefault& t)
+  {
+    vtkm::Particle p = this->Particle.Get(idx);
+    p.Pos = pt;
+    p.NumSteps++;
+    p.Time += t;
+    //std::cout<<idx<<" take step: "<<pt<<" "<<p.NumSteps<<" max= "<<this->MaxSteps<<std::endl;
+    if (p.NumSteps == this->MaxSteps)
+    {
+      p.Status &= ~static_cast<vtkm::Id>(ParticleStatus::SUCCESS);
+      p.Status |= static_cast<vtkm::Id>(ParticleStatus::TERMINATED);
+    }
+    this->Particle.Set(idx, p);
+  }
+
+  /* Set/Change Status */
+  VTKM_EXEC
+  void SetOK(const vtkm::Id& idx)
+  {
+    Clear(idx);
+    SetBit(idx, ParticleStatus::SUCCESS);
+  }
+  VTKM_EXEC
+  void SetTerminated(const vtkm::Id& idx)
+  {
+    ClearBit(idx, ParticleStatus::SUCCESS);
+    SetBit(idx, ParticleStatus::TERMINATED);
+  }
+  VTKM_EXEC
+  void SetTookAnySteps(const vtkm::Id& idx, const bool& val)
+  {
+    if (val)
+      SetBit(idx, ParticleStatus::TOOK_ANY_STEPS);
+    else
+      ClearBit(idx, ParticleStatus::TOOK_ANY_STEPS);
+  }
+  VTKM_EXEC
+  void SetExitSpatialBoundary(const vtkm::Id& idx)
+  {
+    ClearBit(idx, ParticleStatus::SUCCESS);
+    SetBit(idx, ParticleStatus::EXIT_SPATIAL_BOUNDARY);
+  }
+  VTKM_EXEC
+  void SetExitTemporalBoundary(const vtkm::Id& idx)
+  {
+    ClearBit(idx, ParticleStatus::SUCCESS);
+    SetBit(idx, ParticleStatus::EXIT_TEMPORAL_BOUNDARY);
+  }
+  VTKM_EXEC
+  void SetError(const vtkm::Id& idx)
+  {
+    ClearBit(idx, ParticleStatus::SUCCESS);
+    SetBit(idx, ParticleStatus::FAIL);
+  }
+
+  /* Check Status */
+  VTKM_EXEC
+  bool OK(const vtkm::Id& idx) { return CheckBit(idx, ParticleStatus::SUCCESS); }
+  VTKM_EXEC
+  bool Terminated(const vtkm::Id& idx) { return CheckBit(idx, ParticleStatus::TERMINATED); }
+  VTKM_EXEC
+  bool ExitSpatialBoundary(const vtkm::Id& idx)
+  {
+    return CheckBit(idx, ParticleStatus::EXIT_SPATIAL_BOUNDARY);
+  }
+  VTKM_EXEC
+  bool ExitTemporalBoundary(const vtkm::Id& idx)
+  {
+    return CheckBit(idx, ParticleStatus::EXIT_TEMPORAL_BOUNDARY);
+  }
+  VTKM_EXEC
+  bool Error(const vtkm::Id& idx) { return CheckBit(idx, ParticleStatus::FAIL); }
+  VTKM_EXEC
+  bool Done(const vtkm::Id& idx) { return !Integrateable(idx); }
+  VTKM_EXEC
+  bool Integrateable(const vtkm::Id& idx)
+  {
+    return OK(idx) && !(Terminated(idx) || ExitSpatialBoundary(idx) || ExitTemporalBoundary(idx));
+  }
+
+  /* Bit Operations */
+  VTKM_EXEC
+  void Clear(const vtkm::Id& idx)
+  {
+    vtkm::Particle p = Particle.Get(idx);
+    p.Status = 0;
+    Particle.Set(idx, p);
+  }
+  VTKM_EXEC
+  void SetBit(const vtkm::Id& idx, const ParticleStatus& b)
+  {
+    vtkm::Particle p = Particle.Get(idx);
+    p.Status |= static_cast<vtkm::Id>(b);
+    Particle.Set(idx, p);
+  }
+  VTKM_EXEC
+  void ClearBit(const vtkm::Id& idx, const ParticleStatus& b)
+  {
+    vtkm::Particle p = Particle.Get(idx);
+    p.Status &= ~static_cast<vtkm::Id>(b);
+    Particle.Set(idx, p);
+  }
+  VTKM_EXEC
+  bool CheckBit(const vtkm::Id& idx, const ParticleStatus& b) const
+  {
+    return (Particle.Get(idx).Status & static_cast<vtkm::Id>(b)) != 0;
+  }
+
+  VTKM_EXEC
+  vtkm::Vec3f GetPos(const vtkm::Id& idx) const { return Particle.Get(idx).Pos; }
+  VTKM_EXEC
+  vtkm::Id GetStep(const vtkm::Id& idx) const { return Particle.Get(idx).NumSteps; }
+  VTKM_EXEC
+  vtkm::Id GetStatus(const vtkm::Id& idx) const { return Particle.Get(idx).Status; }
+  VTKM_EXEC
+  vtkm::FloatDefault GetTime(const vtkm::Id& idx) const { return Particle.Get(idx).Time; }
+  //  VTKM_EXEC
+  //  void SetTime(const vtkm::Id& idx, vtkm::FloatDefault time) const { Time.Set(idx, time); }
+
+protected:
+  using ParticlePortal =
+    typename vtkm::cont::ArrayHandle<vtkm::Particle>::template ExecutionTypes<Device>::Portal;
+
+  ParticlePortal Particle;
+  vtkm::Id MaxSteps;
+};
+
+
+class ParticlesAOS : public vtkm::cont::ExecutionObjectBase
+{
+public:
+  template <typename Device>
+  VTKM_CONT vtkm::worklet::particleadvection::ParticleExecutionObjectAOS<Device>
+    PrepareForExecution(Device) const
+  {
+
+    return vtkm::worklet::particleadvection::ParticleExecutionObjectAOS<Device>(this->ParticleArray,
+                                                                                this->MaxSteps);
+  }
+
+  VTKM_CONT
+  ParticlesAOS(vtkm::cont::ArrayHandle<vtkm::Particle>& pArray, const vtkm::Id& maxSteps)
+    : ParticleArray(pArray)
+    , MaxSteps(maxSteps)
+  {
+  }
+
+  ParticlesAOS() {}
+
+protected:
+  vtkm::cont::ArrayHandle<vtkm::Particle> ParticleArray;
+  vtkm::Id MaxSteps;
+};
+
+#if 0
+
+template <typename Device>
+class StateRecordingParticleExecutionObject : public ParticleExecutionObject<Device>
+{
+public:
+  VTKM_EXEC_CONT
+  StateRecordingParticleExecutionObject()
+    : ParticleExecutionObject<Device>()
+    , History()
+    , Length(0)
+    , ValidPoint()
+  {
+  }
+
+  StateRecordingParticleExecutionObject(vtkm::cont::ArrayHandle<vtkm::Vec3f> posArray,
+                                        vtkm::cont::ArrayHandle<vtkm::Vec3f> historyArray,
+                                        vtkm::cont::ArrayHandle<vtkm::Id> stepsArray,
+                                        vtkm::cont::ArrayHandle<vtkm::Id> statusArray,
+                                        vtkm::cont::ArrayHandle<vtkm::FloatDefault> timeArray,
+                                        vtkm::cont::ArrayHandle<vtkm::Id> validPointArray,
+                                        vtkm::Id maxSteps)
+    : ParticleExecutionObject<Device>(posArray, stepsArray, statusArray, timeArray, maxSteps)
+  {
+    Length = maxSteps;
+    vtkm::Id numPos = posArray.GetNumberOfValues();
+    History = historyArray.PrepareForOutput(numPos * Length, Device());
+    ValidPoint = validPointArray.PrepareForInPlace(Device());
+  }
+
+  VTKM_EXEC_CONT
+  void TakeStep(const vtkm::Id& idx, const vtkm::Vec3f& pt)
+  {
+    this->ParticleExecutionObject<Device>::TakeStep(idx, pt);
+
+    //TakeStep incremented the step, so we want the PREV step value.
+    vtkm::Id nSteps = this->Steps.Get(idx) - 1;
+
+    // Update the step for streamline storing portals.
+    // This includes updating the history and the valid points.
+    vtkm::Id loc = idx * Length + nSteps;
+    this->History.Set(loc, pt);
+    this->ValidPoint.Set(loc, 1);
+  }
+
+  vtkm::Vec3f GetHistory(const vtkm::Id& idx, const vtkm::Id& step) const
+  {
+    return this->History.Get(idx * this->Length + step);
+  }
+
+protected:
+  using IdPortal =
+    typename vtkm::cont::ArrayHandle<vtkm::Id>::template ExecutionTypes<Device>::Portal;
+  using PositionPortal =
+    typename vtkm::cont::ArrayHandle<vtkm::Vec3f>::template ExecutionTypes<Device>::Portal;
+
+  PositionPortal History;
+  vtkm::Id Length;
+  IdPortal ValidPoint;
+};
+
+class StateRecordingParticles : vtkm::cont::ExecutionObjectBase
+{
+public:
+  template <typename Device>
+  VTKM_CONT vtkm::worklet::particleadvection::StateRecordingParticleExecutionObject<Device>
+    PrepareForExecution(Device) const
+  {
+    return vtkm::worklet::particleadvection::StateRecordingParticleExecutionObject<Device>(
+      PosArray, HistoryArray, StepsArray, StatusArray, TimeArray, ValidPointArray, MaxSteps);
+  }
+
+  VTKM_CONT
+  StateRecordingParticles(vtkm::cont::ArrayHandle<vtkm::Vec3f>& posArray,
+                          vtkm::cont::ArrayHandle<vtkm::Vec3f>& historyArray,
+                          vtkm::cont::ArrayHandle<vtkm::Id>& stepsArray,
+                          vtkm::cont::ArrayHandle<vtkm::Id>& statusArray,
+                          vtkm::cont::ArrayHandle<vtkm::FloatDefault>& timeArray,
+                          vtkm::cont::ArrayHandle<vtkm::Id>& validPointArray,
+                          const vtkm::Id& maxSteps)
+  {
+    PosArray = posArray;
+    HistoryArray = historyArray;
+    StepsArray = stepsArray;
+    StatusArray = statusArray;
+    TimeArray = timeArray;
+    ValidPointArray = validPointArray;
+    MaxSteps = maxSteps;
+  }
+
+protected:
+  vtkm::cont::ArrayHandle<vtkm::Id> StepsArray;
+  vtkm::cont::ArrayHandle<vtkm::Id> StatusArray;
+  vtkm::cont::ArrayHandle<vtkm::FloatDefault> TimeArray;
+  vtkm::cont::ArrayHandle<vtkm::Id> ValidPointArray;
+  vtkm::cont::ArrayHandle<vtkm::Vec3f> HistoryArray;
+  vtkm::cont::ArrayHandle<vtkm::Vec3f> PosArray;
+  vtkm::Id MaxSteps;
+};
+#endif
+
+} //namespace particleadvection
+} //namespace worklet
+} //namespace vtkm
+
+#endif // vtk_m_worklet_particleadvection_Particles_AOS_h

vtkm/worklet/testing/UnitTestParticleAdvection.cxx

@@ -23,6 +23,10 @@
 
 #include <vtkm/io/writer/VTKDataSetWriter.h>
 
+//timers
+#include <chrono>
+#include <ctime>
+
 namespace
 {
 vtkm::FloatDefault vecData[125 * 3] = {
@@ -646,19 +650,95 @@ void TestParticleStatus()
   auto res = pa.Run(rk4, seedsArray, maxSteps);
   auto statusPortal = res.status.GetPortalConstControl();
 
-  vtkm::Id tookStep0 = statusPortal.Get(0) &
-    static_cast<vtkm::Id>(vtkm::worklet::particleadvection::ParticleStatus::TOOK_ANY_STEPS);
-  vtkm::Id tookStep1 = statusPortal.Get(1) &
-    static_cast<vtkm::Id>(vtkm::worklet::particleadvection::ParticleStatus::TOOK_ANY_STEPS);
+  vtkm::Id tookStep0 =
+    statusPortal.Get(0) & static_cast<vtkm::Id>(vtkm::ParticleStatus::TOOK_ANY_STEPS);
+  vtkm::Id tookStep1 =
+    statusPortal.Get(1) & static_cast<vtkm::Id>(vtkm::ParticleStatus::TOOK_ANY_STEPS);
   VTKM_TEST_ASSERT(tookStep0 != 0, "Particle failed to take any steps");
   VTKM_TEST_ASSERT(tookStep1 == 0, "Particle took a step when it should not have.");
 }
 
+
+double TestParticleAOS(bool doAOS)
+{
+  vtkm::Bounds bounds(0, 1, 0, 1, 0, 1);
+  const vtkm::Id3 dims(5, 5, 5);
+  vtkm::cont::DataSet ds = CreateUniformDataSet(bounds, dims);
+
+  vtkm::Id nElements = dims[0] * dims[1] * dims[2] * 3;
+
+  std::vector<vtkm::Vec<vtkm::FloatDefault, 3>> field;
+  for (vtkm::Id i = 0; i < nElements; i++)
+  {
+    vtkm::FloatDefault x = vecData[i];
+    vtkm::FloatDefault y = vecData[++i];
+    vtkm::FloatDefault z = vecData[++i];
+    vtkm::Vec3f vec(x, y, z);
+    field.push_back(vtkm::Normal(vec));
+  }
+
+  vtkm::cont::ArrayHandle<vtkm::Vec3f> fieldArray;
+  fieldArray = vtkm::cont::make_ArrayHandle(field);
+
+  using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
+  using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldHandle>;
+  using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
+  vtkm::Id maxSteps = 1000;
+  vtkm::FloatDefault stepSize = 0.001f;
+
+  GridEvalType eval(ds.GetCoordinateSystem(), ds.GetCellSet(), fieldArray);
+  RK4Type rk4(eval, stepSize);
+
+  vtkm::worklet::ParticleAdvection pa;
+  std::vector<vtkm::Particle> particles;
+  std::vector<vtkm::Vec3f> pts;
+
+  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);
+  }
+
+  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;
+  }
+  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;
+  }
+
+  return dT.count();
+}
+
 void TestParticleAdvection()
 {
-  TestEvaluators();
-  TestParticleWorklets();
-  TestParticleStatus();
+  //  TestEvaluators();
+  //  TestParticleWorklets();
+  //  TestParticleStatus();
+
+  double aosTime = TestParticleAOS(true);
+  double soaTime = TestParticleAOS(false);
+
+  std::cout << "AOS_time= " << aosTime << std::endl;
+  std::cout << "SOA_time= " << soaTime << std::endl;
 }
 
 int UnitTestParticleAdvection(int argc, char* argv[])