warpx particle filter

src/libs/ascent/runtimes/flow_filters/ascent_runtime_filters.cpp

@@ -140,6 +140,7 @@ register_builtin()
     AscentRuntime::register_filter_type<VTKHTriangulate>("transforms","triangulate");
     AscentRuntime::register_filter_type<VTKHParticleAdvection>("transforms","particle_advection");
     AscentRuntime::register_filter_type<VTKHStreamline>("transforms","streamline");
+    AscentRuntime::register_filter_type<VTKHWarpXStreamline>("transforms","warpx_streamline");
 
     AscentRuntime::register_filter_type<RoverXRay>("extracts", "xray");
     AscentRuntime::register_filter_type<RoverVolume>("extracts", "volume");

src/libs/ascent/runtimes/flow_filters/ascent_runtime_vtkh_filters.cpp

@@ -65,6 +65,7 @@
 #include <vtkh/filters/Slice.hpp>
 #include <vtkh/filters/Statistics.hpp>
 #include <vtkh/filters/Streamline.hpp>
+#include <vtkh/filters/WarpXStreamline.hpp>
 #include <vtkh/filters/Threshold.hpp>
 #include <vtkh/filters/Triangulate.hpp>
 #include <vtkh/filters/VectorMagnitude.hpp>
@@ -3927,6 +3928,148 @@ VTKHStreamline::~VTKHStreamline()
 // empty
 }
 
+//-----------------------------------------------------------------------------
+
+VTKHWarpXStreamline::VTKHWarpXStreamline()
+:Filter()
+{
+// empty
+}
+
+//-----------------------------------------------------------------------------
+VTKHWarpXStreamline::~VTKHWarpXStreamline()
+{
+// empty
+}
+
+//-----------------------------------------------------------------------------
+void
+VTKHWarpXStreamline::declare_interface(Node &i)
+{
+    i["type_name"]   = "vtkh_warpx_streamline";
+    i["port_names"].append() = "in";
+    i["output_port"] = "true";
+}
+
+//-----------------------------------------------------------------------------
+bool
+VTKHWarpXStreamline::verify_params(const conduit::Node &params,
+                                     conduit::Node &info)
+{
+    info.reset();
+    bool res = check_string("b_field", params, info, false);
+    res &= check_string("e_field", params, info, false);
+    res &= check_numeric("num_steps", params, info, true, true);
+    res &= check_numeric("step_size", params, info, true, true);
+
+    std::vector<std::string> valid_paths;
+    valid_paths.push_back("b_field");
+    valid_paths.push_back("e_field");
+    valid_paths.push_back("charge_field");
+    valid_paths.push_back("mass_field");
+    valid_paths.push_back("momentum_field");
+    valid_paths.push_back("weighting_field");
+    valid_paths.push_back("num_steps");
+    valid_paths.push_back("step_size");
+
+    std::string surprises = surprise_check(valid_paths, params);
+
+    if(surprises != "")
+    {
+      res = false;
+      info["errors"].append() = surprises;
+    }
+
+    return res;
+}
+//-----------------------------------------------------------------------------
+void
+VTKHWarpXStreamline::execute()
+{
+
+    if(!input(0).check_type<DataObject>())
+    {
+        ASCENT_ERROR("vtkh_warpx_streamline input must be a data object");
+    }
+
+    // grab the data collection and ask for a vtkh collection
+    // which is one vtkh data set per topology
+    DataObject *data_object = input<DataObject>(0);
+    if(!data_object->is_valid())
+    {
+      set_output<DataObject>(data_object);
+      return;
+    }
+    std::shared_ptr<VTKHCollection> collection = data_object->as_vtkh_collection();
+
+    std::string b_field = "B";
+    std::string e_field = "E";
+    std::string charge_field = "Charge";
+    std::string mass_field = "Mass";
+    std::string momentum_field = "Momentum";
+    std::string weighting_field = "Weighting";
+    if(params().has_path("b_field"))
+      b_field = params()["b_field"].as_string();
+    if(params().has_path("e_field"))
+      e_field = params()["e_field"].as_string();
+    if(params().has_path("charge_field"))
+      charge_field = params()["charge_field"].as_string();
+    if(params().has_path("mass_field"))
+      mass_field = params()["mass_field"].as_string();
+    if(params().has_path("momentum_field"))
+      momentum_field = params()["momentum_field"].as_string();
+    if(params().has_path("weighting_field"))
+      weighting_field = params()["weighting_field"].as_string();
+
+    if(!collection->has_field(b_field))
+    {
+      bool throw_error = false;
+      detail::field_error(b_field, this->name(), collection, throw_error);
+      // this creates a data object with an invalid soource
+      set_output<DataObject>(new DataObject());
+      return;
+    }
+    if(!collection->has_field(e_field))
+    {
+      bool throw_error = false;
+      detail::field_error(e_field, this->name(), collection, throw_error);
+      // this creates a data object with an invalid soource
+      set_output<DataObject>(new DataObject());
+      return;
+    }
+
+    std::string topo_name = collection->field_topology(b_field);
+    vtkh::DataSet &data = collection->dataset_by_topology(topo_name);
+
+
+    int numSteps = get_int32(params()["num_steps"], data_object);
+    float stepSize = get_float32(params()["step_size"], data_object);
+
+
+    vtkh::DataSet *output = nullptr;
+    vtkh::WarpXStreamline sl;
+    sl.SetStepSize(stepSize);
+    sl.SetNumberOfSteps(numSteps);
+    sl.SetBField(b_field);
+    sl.SetEField(e_field);
+    sl.SetChargeField(charge_field);
+    sl.SetMassField(mass_field);
+    sl.SetMomentumField(momentum_field);
+    sl.SetWeightingField(weighting_field);
+    sl.SetInput(&data);
+    sl.Update();
+    output = sl.GetOutput();
+
+    // we need to pass through the rest of the topologies, untouched,
+    // and add the result of this operation
+    VTKHCollection *new_coll = collection->copy_without_topology(topo_name);
+    new_coll->add(*output, topo_name);
+    // re wrap in data object
+    DataObject *res =  new DataObject(new_coll);
+    delete output;
+    set_output<DataObject>(res);
+}
+
 //-----------------------------------------------------------------------------
 };
 //-----------------------------------------------------------------------------

src/libs/ascent/runtimes/flow_filters/ascent_runtime_vtkh_filters.hpp

@@ -446,6 +446,21 @@ class ASCENT_API VTKHStreamline : public VTKHParticleAdvection
   virtual void   declare_interface(conduit::Node &i);
 };
 
+//-----------------------------------------------------------------------------
+class ASCENT_API VTKHWarpXStreamline : public ::flow::Filter
+{
+public:
+    VTKHWarpXStreamline();
+    virtual ~VTKHWarpXStreamline();
+
+    virtual void   declare_interface(conduit::Node &i);
+    virtual bool   verify_params(const conduit::Node &params,
+                                 conduit::Node &info);
+    virtual void   execute();
+
+protected:
+    bool record_trajectories = false;
+};
 };
 //-----------------------------------------------------------------------------
 // -- end ascent::runtime::filters --

src/libs/vtkh/filters/CMakeLists.txt

@@ -32,6 +32,7 @@ set(vtkh_filters_headers
     Slice.hpp
     VectorComponent.hpp
     VectorMagnitude.hpp
+    WarpXStreamline.hpp
     )
 
 if(ENABLE_FILTER_CONTOUR_TREE)
@@ -68,6 +69,7 @@ set(vtkh_filters_sources
     Streamline.cpp
     VectorComponent.cpp
     VectorMagnitude.cpp
+    WarpXStreamline.cpp
     )
 
 if(VTKH_ENABLE_FILTER_CONTOUR_TREE)

src/libs/vtkh/filters/WarpXStreamline.cpp

@@ -0,0 +1,188 @@
+#include <iostream>
+#include <vtkh/filters/WarpXStreamline.hpp>
+#include <vtkm/filter/flow/WarpXStreamline.h>
+#include <vtkm/cont/EnvironmentTracker.h>
+#include <vtkh/vtkh.hpp>
+#include <vtkh/Error.hpp>
+
+#if VTKH_PARALLEL
+#include <vtkm/thirdparty/diy/diy.h>
+#include <vtkm/thirdparty/diy/mpi-cast.h>
+#include <mpi.h>
+#endif
+
+
+namespace vtkh
+{
+
+namespace detail
+{
+
+void GenerateChargedParticles(const vtkm::cont::ArrayHandle<vtkm::Vec3f>& pos,
+                              const vtkm::cont::ArrayHandle<vtkm::Vec3f>& mom,
+                              const vtkm::cont::ArrayHandle<vtkm::Float64>& mass,
+                              const vtkm::cont::ArrayHandle<vtkm::Float64>& charge,
+                              const vtkm::cont::ArrayHandle<vtkm::Float64>& weight,
+                              vtkm::cont::ArrayHandle<vtkm::ChargedParticle>& seeds,
+			      const int id_offset)
+{
+  auto pPortal = pos.ReadPortal();
+  auto uPortal = mom.ReadPortal();
+  auto mPortal = mass.ReadPortal();
+  auto qPortal = charge.ReadPortal();
+  auto wPortal = weight.ReadPortal();
+
+  auto numValues = pos.GetNumberOfValues();
+
+  seeds.Allocate(numValues);
+  auto sPortal = seeds.WritePortal();
+
+  for (vtkm::Id i = 0; i < numValues; i++)
+  {
+    vtkm::ChargedParticle electron(
+      pPortal.Get(i), i, mPortal.Get(i), qPortal.Get(i), wPortal.Get(i), uPortal.Get(i));
+    sPortal.Set(i + id_offset, electron);
+  }
+
+}
+
+
+} //end detail
+
+WarpXStreamline::WarpXStreamline()
+: m_e_field_name("E"), 
+	m_b_field_name("B"), 
+	m_charge_field_name("Charge"), 
+	m_mass_field_name("Mass"),
+	m_momentum_field_name("Momentum"), 
+	m_weighting_field_name("Weighting")
+{
+
+}
+
+WarpXStreamline::~WarpXStreamline()
+{
+
+}
+
+void WarpXStreamline::PreExecute()
+{
+  Filter::PreExecute();
+  Filter::CheckForRequiredField(m_b_field_name);
+  Filter::CheckForRequiredField(m_e_field_name);
+  Filter::CheckForRequiredField(m_charge_field_name);
+  Filter::CheckForRequiredField(m_mass_field_name);
+  Filter::CheckForRequiredField(m_momentum_field_name);
+  Filter::CheckForRequiredField(m_weighting_field_name);
+
+}
+
+void WarpXStreamline::PostExecute()
+{
+  Filter::PostExecute();
+}
+
+void WarpXStreamline::DoExecute()
+{
+  this->m_output = new DataSet();
+
+#ifndef VTKH_BYPASS_VTKM_BIH
+
+#ifdef VTKH_PARALLEL
+  // Setup VTK-h and VTK-m comm.
+  MPI_Comm mpi_comm = MPI_Comm_f2c(vtkh::GetMPICommHandle());
+  vtkm::cont::EnvironmentTracker::SetCommunicator(vtkmdiy::mpi::communicator(vtkmdiy::mpi::make_DIY_MPI_Comm(mpi_comm)));
+#endif
+
+  //Make sure that the E field exists on any domain.
+  if (!this->m_input->GlobalFieldExists(m_e_field_name))
+  {
+    throw Error("Domain does not contain specified E vector field for WarpXStreamline analysis.");
+  }
+  //Make sure that the B field exists on any domain.
+  if (!this->m_input->GlobalFieldExists(m_b_field_name))
+  {
+    throw Error("Domain does not contain specified B vector field for WarpXStreamline analysis.");
+  }
+
+  vtkm::cont::PartitionedDataSet inputs;
+
+  vtkm::cont::ArrayHandle<vtkm::ChargedParticle> seeds;
+  //Create charged particles for all domains with the particle spec fields.
+  //TODO: user specified momentum,mass,charge,weighting?
+  if (this->m_input->FieldExists(m_momentum_field_name))
+  {
+    const int num_domains = this->m_input->GetNumberOfDomains();
+    int id_offset = 0;
+    for (int i = 0; i < num_domains; i++)
+    {
+      vtkm::Id domain_id;
+      vtkm::cont::DataSet dom;
+      this->m_input->GetDomain(i, dom, domain_id);
+      if(dom.HasField(m_momentum_field_name))
+      {
+        vtkm::cont::ArrayHandle<vtkm::Vec3f> pos, mom;
+        vtkm::cont::ArrayHandle<vtkm::Float64> mass, charge, w;
+        dom.GetCoordinateSystem().GetData().AsArrayHandle(pos);
+        dom.GetField(m_momentum_field_name).GetData().AsArrayHandle(mom);
+        dom.GetField(m_mass_field_name).GetData().AsArrayHandle(mass);
+        dom.GetField(m_charge_field_name).GetData().AsArrayHandle(charge);
+        dom.GetField(m_weighting_field_name).GetData().AsArrayHandle(w);
+	detail::GenerateChargedParticles(pos,mom,mass,charge,w,seeds, id_offset);
+	//Actual: local unique ids
+	//Question: do we global unique ids? 
+	id_offset += pos.GetNumberOfValues();
+      }
+      if(dom.HasField(m_b_field_name))
+      {
+        auto field = dom.GetField(m_b_field_name).GetData();
+        inputs.AppendPartition(dom);
+      }
+    }
+  }
+  else
+  {
+    throw Error("Domain is missing one or more neccessary fields to create a charged particle: \
+		   Charge, Mass, Momentum, and/or Weighting.");
+  }
+
+  bool validField = (inputs.GetNumberOfPartitions() > 0);
+//Don't really need this check
+//since we got rid of the other check
+#ifdef VTKH_PARALLEL
+  int localNum = static_cast<int>(inputs.GetNumberOfPartitions());
+  int globalNum = 0;
+  MPI_Allreduce((void *)(&localNum),
+                (void *)(&globalNum),
+                1,
+                MPI_INT,
+                MPI_SUM,
+                mpi_comm);
+  validField = (globalNum > 0);
+#endif
+
+  if (!validField)
+  {
+    throw Error("Vector field type does not match a supportable type.");
+  }
+
+  //Everything is valid. Call the VTKm filter.
+
+  vtkm::filter::flow::WarpXStreamline warpxStreamlineFilter;
+
+  warpxStreamlineFilter.SetStepSize(m_step_size);
+  warpxStreamlineFilter.SetBField(m_b_field_name);
+  warpxStreamlineFilter.SetEField(m_e_field_name);
+  warpxStreamlineFilter.SetSeeds(seeds);
+  warpxStreamlineFilter.SetNumberOfSteps(m_num_steps);
+  auto out = warpxStreamlineFilter.Execute(inputs);
+  //out.PrintSummary(std::cerr);
+  int num_domains = m_output->GetNumberOfDomains();
+  for (vtkm::Id i = 0; i < out.GetNumberOfPartitions(); i++)
+  {
+    this->m_output->AddDomain(out.GetPartition(i), num_domains + i);
+  }
+#endif
+}
+
+} //  namespace vtkh