From 2956713174c126dafb7370e7536b0d4fc3d29f37 Mon Sep 17 00:00:00 2001 From: Dave Pugmire Date: Thu, 11 Nov 2021 09:21:15 -0500 Subject: [PATCH] Filters for charged particle advection. --- vtkm/filter/FilterParticleAdvection.h | 21 ++++++++-- vtkm/filter/FilterParticleAdvection.hxx | 17 +++++---- vtkm/filter/FilterTemporalParticleAdvection.h | 7 ++-- .../FilterTemporalParticleAdvection.hxx | 22 ++++++----- vtkm/filter/ParticleAdvection.h | 3 +- vtkm/filter/ParticleAdvection.hxx | 36 +++++++++++++++--- vtkm/filter/PathParticle.h | 3 +- vtkm/filter/PathParticle.hxx | 2 +- vtkm/filter/Pathline.h | 2 +- vtkm/filter/Pathline.hxx | 2 +- vtkm/filter/Streamline.h | 2 +- vtkm/filter/Streamline.hxx | 2 +- .../particleadvection/DataSetIntegrator.h | 25 ++++++++++++ .../particleadvection/DataSetIntegrator.hxx | 16 ++++++++ .../ParticleAdvectionAlgorithm.h | 15 ++++++++ .../testing/UnitTestStreamlineFilter.cxx | 38 +++++++++++++++++++ 16 files changed, 175 insertions(+), 38 deletions(-) diff --git a/vtkm/filter/FilterParticleAdvection.h b/vtkm/filter/FilterParticleAdvection.h index 6e499135f..fdde85f15 100644 --- a/vtkm/filter/FilterParticleAdvection.h +++ b/vtkm/filter/FilterParticleAdvection.h @@ -25,7 +25,7 @@ namespace filter /// Takes as input a vector field and seed locations and advects the seeds /// through the flow field. -template +template class FilterParticleAdvection : public vtkm::filter::FilterDataSetWithField { public: @@ -34,6 +34,16 @@ public: VTKM_CONT FilterParticleAdvection(); + VTKM_CONT + void SetActiveFields( + const std::string& name1, + const std::string& name2, + vtkm::cont::Field::Association association = vtkm::cont::Field::Association::ANY) + { + this->FilterDataSetWithField::SetActiveField(name1, association); + this->ActiveFieldName2 = name2; + } + VTKM_CONT void SetStepSize(vtkm::FloatDefault s) { this->StepSize = s; } @@ -41,14 +51,14 @@ public: void SetNumberOfSteps(vtkm::Id n) { this->NumberOfSteps = n; } VTKM_CONT - void SetSeeds(const std::vector& seeds, + void SetSeeds(const std::vector& seeds, vtkm::CopyFlag copyFlag = vtkm::CopyFlag::On) { this->Seeds = vtkm::cont::make_ArrayHandle(seeds, copyFlag); } VTKM_CONT - void SetSeeds(vtkm::cont::ArrayHandle& seeds) { this->Seeds = seeds; } + void SetSeeds(vtkm::cont::ArrayHandle& seeds) { this->Seeds = seeds; } VTKM_CONT bool GetUseThreadedAlgorithm() { return this->UseThreadedAlgorithm; } @@ -72,13 +82,16 @@ protected: VTKM_CONT virtual void ValidateOptions() const; using DSIType = vtkm::filter::particleadvection::DataSetIntegrator; + using DSIType2 = vtkm::filter::particleadvection::DataSetIntegrator2; + VTKM_CONT std::vector CreateDataSetIntegrators( const vtkm::cont::PartitionedDataSet& input, const vtkm::filter::particleadvection::BoundsMap& boundsMap) const; + std::string ActiveFieldName2; vtkm::Id NumberOfSteps; vtkm::FloatDefault StepSize; - vtkm::cont::ArrayHandle Seeds; + vtkm::cont::ArrayHandle Seeds; bool UseThreadedAlgorithm; private: diff --git a/vtkm/filter/FilterParticleAdvection.hxx b/vtkm/filter/FilterParticleAdvection.hxx index d4f9bf800..3307d8054 100644 --- a/vtkm/filter/FilterParticleAdvection.hxx +++ b/vtkm/filter/FilterParticleAdvection.hxx @@ -18,8 +18,8 @@ namespace filter { //----------------------------------------------------------------------------- -template -inline VTKM_CONT FilterParticleAdvection::FilterParticleAdvection() +template +inline VTKM_CONT FilterParticleAdvection::FilterParticleAdvection() : vtkm::filter::FilterDataSetWithField() , NumberOfSteps(0) , StepSize(0) @@ -27,8 +27,8 @@ inline VTKM_CONT FilterParticleAdvection::FilterParticleAdvection() { } -template -void FilterParticleAdvection::ValidateOptions() const +template +void FilterParticleAdvection::ValidateOptions() const { if (this->GetUseCoordinateSystemAsField()) throw vtkm::cont::ErrorFilterExecution("Coordinate system as field not supported"); @@ -40,9 +40,9 @@ void FilterParticleAdvection::ValidateOptions() const throw vtkm::cont::ErrorFilterExecution("Step size not specified."); } -template +template std::vector -FilterParticleAdvection::CreateDataSetIntegrators( +FilterParticleAdvection::CreateDataSetIntegrators( const vtkm::cont::PartitionedDataSet& input, const vtkm::filter::particleadvection::BoundsMap& boundsMap) const { @@ -66,9 +66,10 @@ FilterParticleAdvection::CreateDataSetIntegrators( } //----------------------------------------------------------------------------- -template +template template -inline VTKM_CONT vtkm::cont::DataSet FilterParticleAdvection::PrepareForExecution( +inline VTKM_CONT vtkm::cont::DataSet +FilterParticleAdvection::PrepareForExecution( const vtkm::cont::DataSet& input, vtkm::filter::PolicyBase policy) { diff --git a/vtkm/filter/FilterTemporalParticleAdvection.h b/vtkm/filter/FilterTemporalParticleAdvection.h index be14b8acf..a5e97aa5a 100644 --- a/vtkm/filter/FilterTemporalParticleAdvection.h +++ b/vtkm/filter/FilterTemporalParticleAdvection.h @@ -25,8 +25,9 @@ namespace filter /// Takes as input a vector field and seed locations and advects the seeds /// through the flow field. -template -class FilterTemporalParticleAdvection : public vtkm::filter::FilterParticleAdvection +template +class FilterTemporalParticleAdvection + : public vtkm::filter::FilterParticleAdvection { public: VTKM_CONT @@ -52,7 +53,7 @@ public: protected: VTKM_CONT void ValidateOptions(const vtkm::cont::PartitionedDataSet& input) const; - using vtkm::filter::FilterParticleAdvection::ValidateOptions; + using vtkm::filter::FilterParticleAdvection::ValidateOptions; using DSIType = vtkm::filter::particleadvection::TemporalDataSetIntegrator; VTKM_CONT std::vector CreateDataSetIntegrators( diff --git a/vtkm/filter/FilterTemporalParticleAdvection.hxx b/vtkm/filter/FilterTemporalParticleAdvection.hxx index 9d9769ce2..eac02971b 100644 --- a/vtkm/filter/FilterTemporalParticleAdvection.hxx +++ b/vtkm/filter/FilterTemporalParticleAdvection.hxx @@ -18,19 +18,20 @@ namespace filter { //----------------------------------------------------------------------------- -template -inline VTKM_CONT FilterTemporalParticleAdvection::FilterTemporalParticleAdvection() - : vtkm::filter::FilterParticleAdvection() +template +inline VTKM_CONT +FilterTemporalParticleAdvection::FilterTemporalParticleAdvection() + : vtkm::filter::FilterParticleAdvection() , PreviousTime(0) , NextTime(0) { } -template -void FilterTemporalParticleAdvection::ValidateOptions( +template +void FilterTemporalParticleAdvection::ValidateOptions( const vtkm::cont::PartitionedDataSet& input) const { - this->vtkm::filter::FilterParticleAdvection::ValidateOptions(); + this->vtkm::filter::FilterParticleAdvection::ValidateOptions(); if (this->NextDataSet.GetNumberOfPartitions() != input.GetNumberOfPartitions()) throw vtkm::cont::ErrorFilterExecution("Number of partitions do not match"); @@ -38,9 +39,9 @@ void FilterTemporalParticleAdvection::ValidateOptions( throw vtkm::cont::ErrorFilterExecution("Previous time must be less than Next time."); } -template +template std::vector -FilterTemporalParticleAdvection::CreateDataSetIntegrators( +FilterTemporalParticleAdvection::CreateDataSetIntegrators( const vtkm::cont::PartitionedDataSet& input, const vtkm::filter::particleadvection::BoundsMap& boundsMap) const { @@ -65,9 +66,10 @@ FilterTemporalParticleAdvection::CreateDataSetIntegrators( } //----------------------------------------------------------------------------- -template +template template -inline VTKM_CONT vtkm::cont::DataSet FilterTemporalParticleAdvection::PrepareForExecution( +inline VTKM_CONT vtkm::cont::DataSet +FilterTemporalParticleAdvection::PrepareForExecution( const vtkm::cont::DataSet& input, vtkm::filter::PolicyBase policy) { diff --git a/vtkm/filter/ParticleAdvection.h b/vtkm/filter/ParticleAdvection.h index 6f6948c1c..6e0bf1f2d 100644 --- a/vtkm/filter/ParticleAdvection.h +++ b/vtkm/filter/ParticleAdvection.h @@ -23,7 +23,8 @@ namespace filter /// Takes as input a vector field and seed locations and generates the /// end points for each seed through the vector field. -class ParticleAdvection : public vtkm::filter::FilterParticleAdvection +class ParticleAdvection + : public vtkm::filter::FilterParticleAdvection { public: VTKM_CONT ParticleAdvection(); diff --git a/vtkm/filter/ParticleAdvection.hxx b/vtkm/filter/ParticleAdvection.hxx index b895ac2a4..1e769adaf 100644 --- a/vtkm/filter/ParticleAdvection.hxx +++ b/vtkm/filter/ParticleAdvection.hxx @@ -23,7 +23,7 @@ namespace filter //----------------------------------------------------------------------------- inline VTKM_CONT ParticleAdvection::ParticleAdvection() - : vtkm::filter::FilterParticleAdvection() + : vtkm::filter::FilterParticleAdvection() { } @@ -34,18 +34,42 @@ inline VTKM_CONT vtkm::cont::PartitionedDataSet ParticleAdvection::PrepareForExe const vtkm::filter::PolicyBase&) { using AlgorithmType = vtkm::filter::particleadvection::ParticleAdvectionAlgorithm; + using AlgorithmType2 = vtkm::filter::particleadvection::ParticleAdvectionAlgorithm2; using ThreadedAlgorithmType = vtkm::filter::particleadvection::ParticleAdvectionThreadedAlgorithm; this->ValidateOptions(); vtkm::filter::particleadvection::BoundsMap boundsMap(input); - auto dsi = this->CreateDataSetIntegrators(input, boundsMap); - if (this->GetUseThreadedAlgorithm()) - return vtkm::filter::particleadvection::RunAlgo( - boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds); + if (this->ActiveFieldName2.empty()) + { + auto dsi = this->CreateDataSetIntegrators(input, boundsMap); + + if (this->GetUseThreadedAlgorithm()) + return vtkm::filter::particleadvection::RunAlgo( + boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds); + else + return vtkm::filter::particleadvection::RunAlgo( + boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds); + } else - return vtkm::filter::particleadvection::RunAlgo( + { + std::vector dsi; + std::string field1 = this->GetActiveFieldName(); + std::string field2 = this->ActiveFieldName2; + + for (vtkm::Id i = 0; i < input.GetNumberOfPartitions(); i++) + { + vtkm::Id blockId = boundsMap.GetLocalBlockId(i); + auto ds = input.GetPartition(i); + if (!ds.HasPointField(field1) || !ds.HasPointField(field2)) + throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation"); + dsi.push_back(DSIType2(ds, blockId, field1, field2)); + } + + std::cout << "Electromagnietc RunAlgo" << std::endl; + return vtkm::filter::particleadvection::RunAlgo( boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds); + } } } diff --git a/vtkm/filter/PathParticle.h b/vtkm/filter/PathParticle.h index 50180c644..e3276f7c4 100644 --- a/vtkm/filter/PathParticle.h +++ b/vtkm/filter/PathParticle.h @@ -21,7 +21,8 @@ namespace filter /// Takes as input a vector field and seed locations and generates the /// paths taken by the seeds through the vector field. -class PathParticle : public vtkm::filter::FilterTemporalParticleAdvection +class PathParticle + : public vtkm::filter::FilterTemporalParticleAdvection { public: VTKM_CONT diff --git a/vtkm/filter/PathParticle.hxx b/vtkm/filter/PathParticle.hxx index 6760b2775..ba0fd878f 100644 --- a/vtkm/filter/PathParticle.hxx +++ b/vtkm/filter/PathParticle.hxx @@ -24,7 +24,7 @@ namespace filter //----------------------------------------------------------------------------- inline VTKM_CONT PathParticle::PathParticle() - : vtkm::filter::FilterTemporalParticleAdvection() + : vtkm::filter::FilterTemporalParticleAdvection() { } diff --git a/vtkm/filter/Pathline.h b/vtkm/filter/Pathline.h index 3ed72619c..49e88061d 100644 --- a/vtkm/filter/Pathline.h +++ b/vtkm/filter/Pathline.h @@ -21,7 +21,7 @@ namespace filter /// Takes as input a vector field and seed locations and generates the /// paths taken by the seeds through the vector field. -class Pathline : public vtkm::filter::FilterTemporalParticleAdvection +class Pathline : public vtkm::filter::FilterTemporalParticleAdvection { public: VTKM_CONT diff --git a/vtkm/filter/Pathline.hxx b/vtkm/filter/Pathline.hxx index 693f8a9e3..1eb392c13 100644 --- a/vtkm/filter/Pathline.hxx +++ b/vtkm/filter/Pathline.hxx @@ -24,7 +24,7 @@ namespace filter //----------------------------------------------------------------------------- inline VTKM_CONT Pathline::Pathline() - : vtkm::filter::FilterTemporalParticleAdvection() + : vtkm::filter::FilterTemporalParticleAdvection() { } diff --git a/vtkm/filter/Streamline.h b/vtkm/filter/Streamline.h index eb4bbbc55..87da916a3 100644 --- a/vtkm/filter/Streamline.h +++ b/vtkm/filter/Streamline.h @@ -21,7 +21,7 @@ namespace filter /// Takes as input a vector field and seed locations and generates the /// paths taken by the seeds through the vector field. -class Streamline : public vtkm::filter::FilterParticleAdvection +class Streamline : public vtkm::filter::FilterParticleAdvection { public: VTKM_CONT diff --git a/vtkm/filter/Streamline.hxx b/vtkm/filter/Streamline.hxx index d37064330..4db49042a 100644 --- a/vtkm/filter/Streamline.hxx +++ b/vtkm/filter/Streamline.hxx @@ -23,7 +23,7 @@ namespace filter //----------------------------------------------------------------------------- inline VTKM_CONT Streamline::Streamline() - : vtkm::filter::FilterParticleAdvection() + : vtkm::filter::FilterParticleAdvection() { } diff --git a/vtkm/filter/particleadvection/DataSetIntegrator.h b/vtkm/filter/particleadvection/DataSetIntegrator.h index 793c7ed9c..439eff6ce 100644 --- a/vtkm/filter/particleadvection/DataSetIntegrator.h +++ b/vtkm/filter/particleadvection/DataSetIntegrator.h @@ -92,6 +92,7 @@ public: DataSetIntegrator(const vtkm::cont::DataSet& ds, vtkm::Id id, const std::string& fieldNm) : DataSetIntegratorBase>>(false, id) + // : DataSetIntegratorBase>(false, id) { using Association = vtkm::cont::Field::Association; using FieldType = vtkm::worklet::particleadvection::VelocityField; @@ -103,6 +104,30 @@ public: } }; +class VTKM_ALWAYS_EXPORT DataSetIntegrator2 + : public DataSetIntegratorBase>>> +{ +public: + DataSetIntegrator2(const vtkm::cont::DataSet& ds, + vtkm::Id id, + const std::string& fieldNm1, + const std::string& fieldNm2) + : DataSetIntegratorBase>>(false, id) + { + using Association = vtkm::cont::Field::Association; + using FieldType = vtkm::worklet::particleadvection::ElectroMagneticField; + using EvalType = vtkm::worklet::particleadvection::GridEvaluator; + Association association = ds.GetField(fieldNm1).GetAssociation(); + auto fieldArray1 = this->GetFieldHandle(ds, fieldNm1); + auto fieldArray2 = this->GetFieldHandle(ds, fieldNm2); + + FieldType field(fieldArray1, fieldArray2, association); + this->Eval = std::shared_ptr(new EvalType(ds, field)); + } +}; + class VTKM_ALWAYS_EXPORT TemporalDataSetIntegrator : public DataSetIntegratorBase>>> diff --git a/vtkm/filter/particleadvection/DataSetIntegrator.hxx b/vtkm/filter/particleadvection/DataSetIntegrator.hxx index da3893277..0ba215e3d 100644 --- a/vtkm/filter/particleadvection/DataSetIntegrator.hxx +++ b/vtkm/filter/particleadvection/DataSetIntegrator.hxx @@ -20,6 +20,8 @@ namespace particleadvection using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator< vtkm::worklet::particleadvection::VelocityField>>; +using GridEvalType2 = vtkm::worklet::particleadvection::GridEvaluator< + vtkm::worklet::particleadvection::ElectroMagneticField>>; using TemporalGridEvalType = vtkm::worklet::particleadvection::TemporalGridEvaluator< vtkm::worklet::particleadvection::VelocityField>>; using RK4Type = vtkm::worklet::particleadvection::RK4Integrator; @@ -42,6 +44,20 @@ inline void DataSetIntegratorBase::DoAdvect( result = Worklet.Run(stepper, seeds, maxSteps); } +//----- +// Specialization for ParticleAdvection worklet w/ electromagnetic field +template <> +template <> +inline void DataSetIntegratorBase::DoAdvect( + vtkm::cont::ArrayHandle& seeds, + const Stepper& stepper, + vtkm::Id maxSteps, + vtkm::worklet::ParticleAdvectionResult& result) const +{ + vtkm::worklet::ParticleAdvection Worklet; + result = Worklet.Run(stepper, seeds, maxSteps); +} + //----- // Specialization for Streamline worklet template <> diff --git a/vtkm/filter/particleadvection/ParticleAdvectionAlgorithm.h b/vtkm/filter/particleadvection/ParticleAdvectionAlgorithm.h index d75cb9237..25998ee3c 100644 --- a/vtkm/filter/particleadvection/ParticleAdvectionAlgorithm.h +++ b/vtkm/filter/particleadvection/ParticleAdvectionAlgorithm.h @@ -23,6 +23,7 @@ namespace particleadvection { using DSIType = vtkm::filter::particleadvection::DataSetIntegrator; +using DSIType2 = vtkm::filter::particleadvection::DataSetIntegrator2; using TDSIType = vtkm::filter::particleadvection::TemporalDataSetIntegrator; //------------------------------------------------------------------------------------------- @@ -41,6 +42,20 @@ public: } }; +//Particle advection +class VTKM_ALWAYS_EXPORT ParticleAdvectionAlgorithm2 + : public AdvectorBaseAlgorithm> +{ +public: + ParticleAdvectionAlgorithm2(const vtkm::filter::particleadvection::BoundsMap& bm, + const std::vector& blocks) + : AdvectorBaseAlgorithm>( + bm, + blocks) + { + } +}; + //Threaded particle advection class VTKM_ALWAYS_EXPORT ParticleAdvectionThreadedAlgorithm : public AdvectorBaseThreadedAlgorithm pts; + vtkm::FloatDefault mass = 1.0f, charge = 1.0f, weight = 1.0; + vtkm::Vec3f momentum(.1, .1, .1); + pts.push_back( + vtkm::ChargedParticle(vtkm::Vec3f(.2f, 1.0f, .2f), 0, mass, charge, weight, momentum)); + pts.push_back( + vtkm::ChargedParticle(vtkm::Vec3f(.2f, 2.0f, .2f), 0, mass, charge, weight, momentum)); + pts.push_back( + vtkm::ChargedParticle(vtkm::Vec3f(.2f, 3.0f, .2f), 0, mass, charge, weight, momentum)); + vtkm::cont::ArrayHandle seedArray = + vtkm::cont::make_ArrayHandle(pts, vtkm::CopyFlag::On); + + /* + vtkm::filter::ParticleAdvection filter; + filter.SetStepSize(0.1f); + filter.SetNumberOfSteps(20); + filter.SetSeeds(seedArray); + + filter.SetActiveFields(fieldName, fieldName); + auto output = filter.Execute(ds); + + //Validate the result is correct. + VTKM_TEST_ASSERT(output.GetNumberOfCoordinateSystems() == 1, + "Wrong number of coordinate systems in the output dataset"); + + vtkm::cont::CoordinateSystem coords = output.GetCoordinateSystem(); + VTKM_TEST_ASSERT(coords.GetNumberOfPoints() == 63, "Wrong number of coordinates"); + + vtkm::cont::DynamicCellSet dcells = output.GetCellSet(); + VTKM_TEST_ASSERT(dcells.GetNumberOfCells() == 3, "Wrong number of cells"); + */ + } } void TestPathline()