From 2dc485ac22d8651d0256447214a981fd2119a1e6 Mon Sep 17 00:00:00 2001 From: Abhishek Yenpure Date: Sun, 14 Aug 2022 15:50:36 -0700 Subject: [PATCH 1/4] Adding WarpX code path for steady state streamlines - Adding implementation for advecting charged particles in a electromagnetic field (This was done as a POC for WarpX w/ Axel Huebl from LBNL) --- .../flow/NewFilterParticleAdvection.cxx | 2 +- vtkm/filter/flow/NewFilterParticleAdvection.h | 7 ++ .../NewFilterParticleAdvectionSteadyState.cxx | 17 +-- vtkm/filter/flow/internal/DataSetIntegrator.h | 4 +- .../internal/DataSetIntegratorSteadyState.h | 100 ++++++++++++++---- .../internal/DataSetIntegratorUnsteadyState.h | 4 +- 6 files changed, 102 insertions(+), 32 deletions(-) diff --git a/vtkm/filter/flow/NewFilterParticleAdvection.cxx b/vtkm/filter/flow/NewFilterParticleAdvection.cxx index 3ef938356..7916876d2 100644 --- a/vtkm/filter/flow/NewFilterParticleAdvection.cxx +++ b/vtkm/filter/flow/NewFilterParticleAdvection.cxx @@ -37,7 +37,7 @@ VTKM_CONT void NewFilterParticleAdvection::ValidateOptions() const if (this->Seeds.GetNumberOfValues() == 0) throw vtkm::cont::ErrorFilterExecution("No seeds provided."); if (!this->Seeds.IsBaseComponentType() && - this->Seeds.IsBaseComponentType()) + !this->Seeds.IsBaseComponentType()) throw vtkm::cont::ErrorFilterExecution("Unsupported particle type in seed array."); if (this->NumberOfSteps == 0) throw vtkm::cont::ErrorFilterExecution("Number of steps not specified."); diff --git a/vtkm/filter/flow/NewFilterParticleAdvection.h b/vtkm/filter/flow/NewFilterParticleAdvection.h index ee18f6dbb..85593400b 100644 --- a/vtkm/filter/flow/NewFilterParticleAdvection.h +++ b/vtkm/filter/flow/NewFilterParticleAdvection.h @@ -51,12 +51,19 @@ public: VTKM_CONT void SetSolverRK4() { this->SolverType = vtkm::filter::flow::IntegrationSolverType::RK4_TYPE; } + VTKM_CONT void SetSolverEuler() { this->SolverType = vtkm::filter::flow::IntegrationSolverType::EULER_TYPE; } + VTKM_CONT + void SetVectorFieldType(vtkm::filter::flow::VectorFieldType vecFieldType) + { + this->VecFieldType = vecFieldType; + } + VTKM_CONT bool GetUseThreadedAlgorithm() { return this->UseThreadedAlgorithm; } diff --git a/vtkm/filter/flow/NewFilterParticleAdvectionSteadyState.cxx b/vtkm/filter/flow/NewFilterParticleAdvectionSteadyState.cxx index fbee0fa1c..f2857c7c8 100644 --- a/vtkm/filter/flow/NewFilterParticleAdvectionSteadyState.cxx +++ b/vtkm/filter/flow/NewFilterParticleAdvectionSteadyState.cxx @@ -25,7 +25,7 @@ namespace { VTKM_CONT std::vector CreateDataSetIntegrators(const vtkm::cont::PartitionedDataSet& input, - const std::string& activeField, + const std::pair& activeField, const vtkm::filter::flow::internal::BoundsMap& boundsMap, const vtkm::filter::flow::IntegrationSolverType solverType, const vtkm::filter::flow::VectorFieldType vecFieldType, @@ -38,7 +38,9 @@ CreateDataSetIntegrators(const vtkm::cont::PartitionedDataSet& input, { vtkm::Id blockId = boundsMap.GetLocalBlockId(i); auto ds = input.GetPartition(i); - if (!ds.HasPointField(activeField) && !ds.HasCellField(activeField)) + if (!ds.HasPointField(activeField.first) && !ds.HasCellField(activeField.first)) + throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation"); + if (!ds.HasPointField(activeField.second) && !ds.HasCellField(activeField.second)) throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation"); dsi.emplace_back(ds, blockId, activeField, solverType, vecFieldType, resultType); @@ -46,6 +48,7 @@ CreateDataSetIntegrators(const vtkm::cont::PartitionedDataSet& input, return dsi; } + } //anonymous namespace VTKM_CONT vtkm::cont::PartitionedDataSet NewFilterParticleAdvectionSteadyState::DoExecutePartitions( @@ -54,13 +57,11 @@ VTKM_CONT vtkm::cont::PartitionedDataSet NewFilterParticleAdvectionSteadyState:: using DSIType = vtkm::filter::flow::internal::DataSetIntegratorSteadyState; this->ValidateOptions(); + std::pair field("E", "B"); + vtkm::filter::flow::internal::BoundsMap boundsMap(input); - auto dsi = CreateDataSetIntegrators(input, - this->GetActiveFieldName(), - boundsMap, - this->SolverType, - this->VecFieldType, - this->GetResultType()); + auto dsi = CreateDataSetIntegrators( + input, field, boundsMap, this->SolverType, this->VecFieldType, this->GetResultType()); vtkm::filter::flow::internal::ParticleAdvector pav( boundsMap, dsi, this->UseThreadedAlgorithm, this->GetResultType()); diff --git a/vtkm/filter/flow/internal/DataSetIntegrator.h b/vtkm/filter/flow/internal/DataSetIntegrator.h index b3974cda7..3f2d21881 100644 --- a/vtkm/filter/flow/internal/DataSetIntegrator.h +++ b/vtkm/filter/flow/internal/DataSetIntegrator.h @@ -60,9 +60,11 @@ using DSIHelperInfoType = vtkm::cont::internal::Variant class DataSetIntegrator { -protected: +public: using VelocityFieldNameType = std::string; using ElectroMagneticFieldNameType = std::pair; + +protected: using FieldNameType = vtkm::cont::internal::Variant; diff --git a/vtkm/filter/flow/internal/DataSetIntegratorSteadyState.h b/vtkm/filter/flow/internal/DataSetIntegratorSteadyState.h index 6e4bd5712..e5828d3a6 100644 --- a/vtkm/filter/flow/internal/DataSetIntegratorSteadyState.h +++ b/vtkm/filter/flow/internal/DataSetIntegratorSteadyState.h @@ -68,6 +68,26 @@ protected: throw vtkm::cont::ErrorFilterExecution("Velocity field vector type not available"); } + template + VTKM_CONT void GetElectroMagneticField( + vtkm::worklet::flow::ElectroMagneticField& ebField) const + { + if (this->FieldName.GetIndex() == this->FieldName.GetIndexOf()) + { + const auto& fieldNm = this->FieldName.Get(); + const auto& electric = fieldNm.first; + const auto& magnetic = fieldNm.second; + auto assoc = this->DataSet.GetField(electric).GetAssociation(); + + ArrayType eField, bField; + vtkm::cont::ArrayCopyShallowIfPossible(this->DataSet.GetField(electric).GetData(), eField); + vtkm::cont::ArrayCopyShallowIfPossible(this->DataSet.GetField(magnetic).GetData(), bField); + ebField = vtkm::worklet::flow::ElectroMagneticField(eField, bField, assoc); + } + else + throw vtkm::cont::ErrorFilterExecution("Velocity field vector type not available"); + } + private: vtkm::cont::DataSet DataSet; }; @@ -76,17 +96,23 @@ private: namespace internal { using ArrayType = vtkm::cont::ArrayHandle; + using VelocityFieldType = vtkm::worklet::flow::VelocityField; -using SteadyStateGridEvalType = vtkm::worklet::flow::GridEvaluator; +using VelocityEvalType = vtkm::worklet::flow::GridEvaluator; + +using EBFieldType = vtkm::worklet::flow::ElectroMagneticField; +using EBEvalType = vtkm::worklet::flow::GridEvaluator; -template -class AdvectHelper; +//template +//class AdvectHelper; -template -class AdvectHelper +template +class AdvectHelper // { public: - static void Advect(const VelocityFieldType& velField, + using SteadyStateGridEvalType = vtkm::worklet::flow::GridEvaluator; + + static void Advect(const FieldType& vecField, const vtkm::cont::DataSet& ds, vtkm::cont::ArrayHandle& seedArray, vtkm::FloatDefault stepSize, @@ -99,20 +125,20 @@ public: DoAdvect( - velField, ds, seedArray, stepSize, maxSteps, result); + vecField, ds, seedArray, stepSize, maxSteps, result); } else if (solverType == IntegrationSolverType::EULER_TYPE) { DoAdvect( - velField, ds, seedArray, stepSize, maxSteps, result); + vecField, ds, seedArray, stepSize, maxSteps, result); } else throw vtkm::cont::ErrorFilterExecution("Unsupported Integrator type"); } - static void Advect(const VelocityFieldType& velField, + static void Advect(const FieldType& vecField, const vtkm::cont::DataSet& ds, vtkm::cont::ArrayHandle& seedArray, vtkm::FloatDefault stepSize, @@ -125,14 +151,14 @@ public: DoAdvect( - velField, ds, seedArray, stepSize, maxSteps, result); + vecField, ds, seedArray, stepSize, maxSteps, result); } else if (solverType == IntegrationSolverType::EULER_TYPE) { DoAdvect( - velField, ds, seedArray, stepSize, maxSteps, result); + vecField, ds, seedArray, stepSize, maxSteps, result); } else throw vtkm::cont::ErrorFilterExecution("Unsupported Integrator type"); @@ -143,7 +169,7 @@ public: class ResultType, template class SolverType> - static void DoAdvect(const VelocityFieldType& velField, + static void DoAdvect(const FieldType& vecField, const vtkm::cont::DataSet& ds, vtkm::cont::ArrayHandle& seedArray, vtkm::FloatDefault stepSize, @@ -154,7 +180,7 @@ public: vtkm::worklet::flow::Stepper, SteadyStateGridEvalType>; WorkletType worklet; - SteadyStateGridEvalType eval(ds, velField); + SteadyStateGridEvalType eval(ds, vecField); StepperType stepper(eval, stepSize); result = worklet.Run(stepper, seedArray, maxSteps); } @@ -173,23 +199,23 @@ VTKM_CONT inline void DataSetIntegratorSteadyState::DoAdvect(DSIHelperInfoVecFieldType == VectorFieldType::VELOCITY_FIELD_TYPE) { using FieldType = vtkm::worklet::flow::VelocityField; - FieldType velField; - this->GetVelocityField(velField); + FieldType vecField; + this->GetVelocityField(vecField); - using AHType = internal::AdvectHelper; + using AHType = internal::AdvectHelper; if (this->IsParticleAdvectionResult()) { vtkm::worklet::flow::ParticleAdvectionResult result; AHType::Advect( - velField, this->DataSet, seedArray, stepSize, maxSteps, this->SolverType, result); + vecField, this->DataSet, seedArray, stepSize, maxSteps, this->SolverType, result); this->UpdateResult(result, b); } else if (this->IsStreamlineResult()) { vtkm::worklet::flow::StreamlineResult result; AHType::Advect( - velField, this->DataSet, seedArray, stepSize, maxSteps, this->SolverType, result); + vecField, this->DataSet, seedArray, stepSize, maxSteps, this->SolverType, result); this->UpdateResult(result, b); } else @@ -200,10 +226,42 @@ VTKM_CONT inline void DataSetIntegratorSteadyState::DoAdvect(DSIHelperInfo& vtkmNotUsed(b), - vtkm::FloatDefault vtkmNotUsed(stepSize), - vtkm::Id vtkmNotUsed(maxSteps)) + DSIHelperInfo& b, + vtkm::FloatDefault stepSize, + vtkm::Id maxSteps) { + using ArrayType = vtkm::cont::ArrayHandle; + + auto copyFlag = (this->CopySeedArray ? vtkm::CopyFlag::On : vtkm::CopyFlag::Off); + auto seedArray = vtkm::cont::make_ArrayHandle(b.V, copyFlag); + + if (this->VecFieldType == VectorFieldType::ELECTRO_MAGNETIC_FIELD_TYPE) + { + using FieldType = vtkm::worklet::flow::ElectroMagneticField; + FieldType ebField; + this->GetElectroMagneticField(ebField); + + using AHType = internal::AdvectHelper; + + if (this->IsParticleAdvectionResult()) + { + vtkm::worklet::flow::ParticleAdvectionResult result; + AHType::Advect( + ebField, this->DataSet, seedArray, stepSize, maxSteps, this->SolverType, result); + this->UpdateResult(result, b); + } + else if (this->IsStreamlineResult()) + { + vtkm::worklet::flow::StreamlineResult result; + AHType::Advect( + ebField, this->DataSet, seedArray, stepSize, maxSteps, this->SolverType, result); + this->UpdateResult(result, b); + } + else + throw vtkm::cont::ErrorFilterExecution("Unsupported result type"); + } + else + throw vtkm::cont::ErrorFilterExecution("Unsupported vector field type"); } } diff --git a/vtkm/filter/flow/internal/DataSetIntegratorUnsteadyState.h b/vtkm/filter/flow/internal/DataSetIntegratorUnsteadyState.h index de4a9d540..0d3aac264 100644 --- a/vtkm/filter/flow/internal/DataSetIntegratorUnsteadyState.h +++ b/vtkm/filter/flow/internal/DataSetIntegratorUnsteadyState.h @@ -254,9 +254,11 @@ VTKM_CONT inline void DataSetIntegratorUnsteadyState::DoAdvect( vtkm::FloatDefault vtkmNotUsed(stepSize), vtkm::Id vtkmNotUsed(maxSteps)) { + throw vtkm::cont::ErrorFilterExecution( + "Unsupported operation : charged particles and electromagnetic fielfs currently only supported " + "for steady state"); } - } } } -- GitLab From 606223edd292243e5b10aeaecb03de9c87ad21cf Mon Sep 17 00:00:00 2001 From: Abhishek Yenpure Date: Sun, 14 Aug 2022 18:58:51 -0700 Subject: [PATCH 2/4] Fixing choice of dataset integrator Adding updates Adding WarpX data and unit test Fixing code from feedback Fixing code from Ollie and Dave's feedback Reducing WarpX dataset size Fixing high precision requirement to store properties Fixing Particle Sizeof Fixing high precision requirement to store properties fixing code from Ollie and Dave's feedback Trying test Fixing ChargedParticle serailization for MPI --- data/data/misc/warpXfields.vtk | 3 + data/data/misc/warpXparticles.vtk | 3 + vtkm/Particle.h | 37 +++-- vtkm/filter/flow/NewFilterParticleAdvection.h | 11 ++ .../NewFilterParticleAdvectionSteadyState.cxx | 54 ++++++-- .../internal/DataSetIntegratorSteadyState.h | 15 +- vtkm/filter/flow/testing/CMakeLists.txt | 1 + .../testing/UnitTestStreamlineFilterWarpX.cxx | 131 ++++++++++++++++++ 8 files changed, 217 insertions(+), 38 deletions(-) create mode 100644 data/data/misc/warpXfields.vtk create mode 100644 data/data/misc/warpXparticles.vtk create mode 100644 vtkm/filter/flow/testing/UnitTestStreamlineFilterWarpX.cxx diff --git a/data/data/misc/warpXfields.vtk b/data/data/misc/warpXfields.vtk new file mode 100644 index 000000000..90bbbda86 --- /dev/null +++ b/data/data/misc/warpXfields.vtk @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:7a569aff0c6872611ef75a4dcb597e1320cb966b80f840a0dbd76a29b2760dbb +size 50335052 diff --git a/data/data/misc/warpXparticles.vtk b/data/data/misc/warpXparticles.vtk new file mode 100644 index 000000000..14b0ff0d8 --- /dev/null +++ b/data/data/misc/warpXparticles.vtk @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:671345bdb045aeadc8e9fa1060de51e53286c74929c6c8c60a529b318f02bbfc +size 3795 diff --git a/vtkm/Particle.h b/vtkm/Particle.h index 027df2878..fc253a14a 100644 --- a/vtkm/Particle.h +++ b/vtkm/Particle.h @@ -164,9 +164,9 @@ public: VTKM_EXEC_CONT ChargedParticle(const vtkm::Vec3f& position, const vtkm::Id& id, - const vtkm::FloatDefault& mass, - const vtkm::FloatDefault& charge, - const vtkm::FloatDefault& weighting, + const vtkm::Float64& mass, + const vtkm::Float64& charge, + const vtkm::Float64& weighting, const vtkm::Vec3f& momentum, const vtkm::Id& numSteps = 0, const vtkm::ParticleStatus& status = vtkm::ParticleStatus(), @@ -184,16 +184,16 @@ public: } VTKM_EXEC_CONT - vtkm::FloatDefault Gamma(vtkm::Vec3f momentum, bool reciprocal = false) const + vtkm::Float64 Gamma(vtkm::Vec3f momentum, bool reciprocal = false) const { constexpr vtkm::FloatDefault c2 = SPEED_OF_LIGHT * SPEED_OF_LIGHT; - const auto fMom2 = vtkm::MagnitudeSquared(momentum); - const auto m2 = this->Mass * this->Mass; - const auto m2_c2_reci = 1.0 / (m2 * c2); + const vtkm::Float64 fMom2 = vtkm::MagnitudeSquared(momentum); + const vtkm::Float64 m2 = this->Mass * this->Mass; + const vtkm::Float64 m2_c2_reci = 1.0 / (m2 * c2); if (reciprocal) - return static_cast(vtkm::RSqrt(1.0 + fMom2 * m2_c2_reci)); + return vtkm::RSqrt(1.0 + fMom2 * m2_c2_reci); else - return static_cast(vtkm::Sqrt(1.0 + fMom2 * m2_c2_reci)); + return vtkm::Sqrt(1.0 + fMom2 * m2_c2_reci); } VTKM_EXEC_CONT @@ -208,11 +208,11 @@ public: vtkm::Vec3f eField = vectors[0]; vtkm::Vec3f bField = vectors[1]; - const vtkm::FloatDefault QoM = this->Charge / this->Mass; + const vtkm::Float64 QoM = this->Charge / this->Mass; const vtkm::Vec3f mom_minus = this->Momentum + (0.5 * this->Charge * eField * length); // Get reciprocal of Gamma - vtkm::Vec3f gamma_reci = this->Gamma(mom_minus, true); + vtkm::Vec3f gamma_reci = static_cast(this->Gamma(mom_minus, true)); const vtkm::Vec3f t = 0.5 * QoM * length * bField * gamma_reci; const vtkm::Vec3f s = 2.0f * t * (1.0 / (1.0 + vtkm::Magnitude(t))); const vtkm::Vec3f mom_prime = mom_minus + vtkm::Cross(mom_minus, t); @@ -254,9 +254,9 @@ public: vtkm::FloatDefault Time = 0; private: - vtkm::FloatDefault Mass; - vtkm::FloatDefault Charge; - vtkm::FloatDefault Weighting; + vtkm::Float64 Mass; + vtkm::Float64 Charge; + vtkm::Float64 Weighting; vtkm::Vec3f Momentum; constexpr static vtkm::FloatDefault SPEED_OF_LIGHT = static_cast(2.99792458e8); @@ -271,11 +271,10 @@ public: + sizeof(vtkm::Id) // NumSteps + sizeof(vtkm::UInt8) // Status + sizeof(vtkm::FloatDefault) // Time - + sizeof(vtkm::FloatDefault) //Mass - + sizeof(vtkm::FloatDefault) //Charge - + sizeof(vtkm::FloatDefault) //Weighting - + sizeof(vtkm::Vec3f) //Momentum - + sizeof(vtkm::FloatDefault); //Speed_of_light + + sizeof(vtkm::Float64) //Mass + + sizeof(vtkm::Float64) //Charge + + sizeof(vtkm::Float64) //Weighting + + sizeof(vtkm::Vec3f); //Momentum return sz; } diff --git a/vtkm/filter/flow/NewFilterParticleAdvection.h b/vtkm/filter/flow/NewFilterParticleAdvection.h index 85593400b..e4bd5ad35 100644 --- a/vtkm/filter/flow/NewFilterParticleAdvection.h +++ b/vtkm/filter/flow/NewFilterParticleAdvection.h @@ -64,6 +64,17 @@ public: this->VecFieldType = vecFieldType; } + //@{ + /// Choose the field to operate on. Note, if + /// `this->UseCoordinateSystemAsField` is true, then the active field is not used. + VTKM_CONT void SetEField(const std::string& name) { this->SetActiveField(0, name); } + + VTKM_CONT void SetBField(const std::string& name) { this->SetActiveField(1, name); } + + VTKM_CONT std::string GetEField() const { return this->GetActiveFieldName(0); } + + VTKM_CONT std::string GetBField() const { return this->GetActiveFieldName(1); } + VTKM_CONT bool GetUseThreadedAlgorithm() { return this->UseThreadedAlgorithm; } diff --git a/vtkm/filter/flow/NewFilterParticleAdvectionSteadyState.cxx b/vtkm/filter/flow/NewFilterParticleAdvectionSteadyState.cxx index f2857c7c8..70c6c27e9 100644 --- a/vtkm/filter/flow/NewFilterParticleAdvectionSteadyState.cxx +++ b/vtkm/filter/flow/NewFilterParticleAdvectionSteadyState.cxx @@ -24,12 +24,14 @@ namespace flow namespace { VTKM_CONT std::vector -CreateDataSetIntegrators(const vtkm::cont::PartitionedDataSet& input, - const std::pair& activeField, - const vtkm::filter::flow::internal::BoundsMap& boundsMap, - const vtkm::filter::flow::IntegrationSolverType solverType, - const vtkm::filter::flow::VectorFieldType vecFieldType, - const vtkm::filter::flow::FlowResultType resultType) +CreateDataSetIntegrators( + const vtkm::cont::PartitionedDataSet& input, + const vtkm::cont::internal::Variant>& + activeField, + const vtkm::filter::flow::internal::BoundsMap& boundsMap, + const vtkm::filter::flow::IntegrationSolverType solverType, + const vtkm::filter::flow::VectorFieldType vecFieldType, + const vtkm::filter::flow::FlowResultType resultType) { using DSIType = vtkm::filter::flow::internal::DataSetIntegratorSteadyState; @@ -38,14 +40,24 @@ CreateDataSetIntegrators(const vtkm::cont::PartitionedDataSet& input, { vtkm::Id blockId = boundsMap.GetLocalBlockId(i); auto ds = input.GetPartition(i); - if (!ds.HasPointField(activeField.first) && !ds.HasCellField(activeField.first)) - throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation"); - if (!ds.HasPointField(activeField.second) && !ds.HasCellField(activeField.second)) - throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation"); - + if (activeField.IsType()) + { + const auto& fieldNm = activeField.Get(); + if (!ds.HasPointField(fieldNm) && !ds.HasCellField(fieldNm)) + throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation"); + } + else if (activeField.IsType()) + { + const auto& fieldNm = activeField.Get(); + const auto& electric = fieldNm.first; + const auto& magnetic = fieldNm.second; + if (!ds.HasPointField(electric) && !ds.HasCellField(electric)) + throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation"); + if (!ds.HasPointField(magnetic) && !ds.HasCellField(magnetic)) + throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation"); + } dsi.emplace_back(ds, blockId, activeField, solverType, vecFieldType, resultType); } - return dsi; } @@ -57,11 +69,25 @@ VTKM_CONT vtkm::cont::PartitionedDataSet NewFilterParticleAdvectionSteadyState:: using DSIType = vtkm::filter::flow::internal::DataSetIntegratorSteadyState; this->ValidateOptions(); - std::pair field("E", "B"); + using VariantType = + vtkm::cont::internal::Variant>; + VariantType variant; + + if (this->VecFieldType == vtkm::filter::flow::VectorFieldType::VELOCITY_FIELD_TYPE) + { + const auto& field = this->GetActiveFieldName(); + variant.Emplace(field); + } + else if (this->VecFieldType == vtkm::filter::flow::VectorFieldType::ELECTRO_MAGNETIC_FIELD_TYPE) + { + const auto& electric = this->GetEField(); + const auto& magnetic = this->GetBField(); + variant.Emplace(electric, magnetic); + } vtkm::filter::flow::internal::BoundsMap boundsMap(input); auto dsi = CreateDataSetIntegrators( - input, field, boundsMap, this->SolverType, this->VecFieldType, this->GetResultType()); + input, variant, boundsMap, this->SolverType, this->VecFieldType, this->GetResultType()); vtkm::filter::flow::internal::ParticleAdvector pav( boundsMap, dsi, this->UseThreadedAlgorithm, this->GetResultType()); diff --git a/vtkm/filter/flow/internal/DataSetIntegratorSteadyState.h b/vtkm/filter/flow/internal/DataSetIntegratorSteadyState.h index e5828d3a6..bd267a044 100644 --- a/vtkm/filter/flow/internal/DataSetIntegratorSteadyState.h +++ b/vtkm/filter/flow/internal/DataSetIntegratorSteadyState.h @@ -55,7 +55,7 @@ protected: VTKM_CONT void GetVelocityField( vtkm::worklet::flow::VelocityField& velocityField) const { - if (this->FieldName.GetIndex() == this->FieldName.GetIndexOf()) + if (this->FieldName.IsType()) { const auto& fieldNm = this->FieldName.Get(); auto assoc = this->DataSet.GetField(fieldNm).GetAssociation(); @@ -72,20 +72,25 @@ protected: VTKM_CONT void GetElectroMagneticField( vtkm::worklet::flow::ElectroMagneticField& ebField) const { - if (this->FieldName.GetIndex() == this->FieldName.GetIndexOf()) + if (this->FieldName.IsType()) { const auto& fieldNm = this->FieldName.Get(); const auto& electric = fieldNm.first; const auto& magnetic = fieldNm.second; - auto assoc = this->DataSet.GetField(electric).GetAssociation(); + auto eAssoc = this->DataSet.GetField(electric).GetAssociation(); + auto bAssoc = this->DataSet.GetField(magnetic).GetAssociation(); + if (eAssoc != bAssoc) + { + throw vtkm::cont::ErrorFilterExecution("E and B field need to have same association"); + } ArrayType eField, bField; vtkm::cont::ArrayCopyShallowIfPossible(this->DataSet.GetField(electric).GetData(), eField); vtkm::cont::ArrayCopyShallowIfPossible(this->DataSet.GetField(magnetic).GetData(), bField); - ebField = vtkm::worklet::flow::ElectroMagneticField(eField, bField, assoc); + ebField = vtkm::worklet::flow::ElectroMagneticField(eField, bField, eAssoc); } else - throw vtkm::cont::ErrorFilterExecution("Velocity field vector type not available"); + throw vtkm::cont::ErrorFilterExecution("Electromagnetic field vector type not available"); } private: diff --git a/vtkm/filter/flow/testing/CMakeLists.txt b/vtkm/filter/flow/testing/CMakeLists.txt index 858aed60a..0e16137fe 100644 --- a/vtkm/filter/flow/testing/CMakeLists.txt +++ b/vtkm/filter/flow/testing/CMakeLists.txt @@ -10,6 +10,7 @@ set(filter_unit_tests UnitTestStreamlineFilter.cxx + UnitTestStreamlineFilterWarpX.cxx UnitTestStreamSurfaceFilter.cxx ) set(worklet_unit_tests diff --git a/vtkm/filter/flow/testing/UnitTestStreamlineFilterWarpX.cxx b/vtkm/filter/flow/testing/UnitTestStreamlineFilterWarpX.cxx new file mode 100644 index 000000000..7813cb90b --- /dev/null +++ b/vtkm/filter/flow/testing/UnitTestStreamlineFilterWarpX.cxx @@ -0,0 +1,131 @@ +//============================================================================ +// 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. +//============================================================================ + +#include + +#include +#include +#include +#include +#include + +namespace +{ + +class GetChargedParticles : public vtkm::worklet::WorkletMapField +{ +public: + GetChargedParticles() {} + + using ControlSignature = void(FieldIn pos, + FieldIn mom, + FieldIn mass, + FieldIn charge, + FieldIn weighting, + FieldOut electrons); + + using ExecutionSignature = void(WorkIndex, _1, _2, _3, _4, _5, _6); + + void operator()(const vtkm::Id index, + const vtkm::Vec3f& pos, + const vtkm::Vec3f& mom, + const vtkm::Float64& mass, + const vtkm::Float64& charge, + const vtkm::Float64& w, + vtkm::ChargedParticle& electron) const + { + electron = vtkm::ChargedParticle(pos, index, mass, charge, w, mom); + } +}; + + +void GenerateChargedParticles(const vtkm::cont::ArrayHandle& pos, + const vtkm::cont::ArrayHandle& mom, + const vtkm::cont::ArrayHandle& mass, + const vtkm::cont::ArrayHandle& charge, + const vtkm::cont::ArrayHandle& weight, + vtkm::cont::ArrayHandle& seeds) +{ + vtkm::cont::Invoker invoker; + invoker(GetChargedParticles{}, pos, mom, mass, charge, weight, seeds); +} + + +void TestStreamlineFilters() +{ + std::string particleFile = vtkm::cont::testing::Testing::DataPath("misc/warpXparticles.vtk"); + std::string fieldFile = vtkm::cont::testing::Testing::DataPath("misc/warpXfields.vtk"); + + using SeedsType = vtkm::cont::ArrayHandle; + + SeedsType seeds; + vtkm::io::VTKDataSetReader seedsReader(particleFile); + vtkm::cont::DataSet seedsData = seedsReader.ReadDataSet(); + vtkm::cont::ArrayHandle pos, mom; + vtkm::cont::ArrayHandle mass, charge, w; + + seedsData.GetCoordinateSystem().GetDataAsDefaultFloat().AsArrayHandle(pos); + seedsData.GetField("Momentum").GetDataAsDefaultFloat().AsArrayHandle(mom); + seedsData.GetField("Mass").GetData().AsArrayHandle(mass); + seedsData.GetField("Charge").GetData().AsArrayHandle(charge); + seedsData.GetField("Weighting").GetData().AsArrayHandle(w); + + GenerateChargedParticles(pos, mom, mass, charge, w, seeds); + + vtkm::io::VTKDataSetReader dataReader(fieldFile); + vtkm::cont::DataSet dataset = dataReader.ReadDataSet(); + vtkm::cont::UnknownCellSet cells = dataset.GetCellSet(); + vtkm::cont::CoordinateSystem coords = dataset.GetCoordinateSystem(); + + auto bounds = coords.GetBounds(); + std::cout << "Bounds : " << bounds << std::endl; + using Structured3DType = vtkm::cont::CellSetStructured<3>; + Structured3DType castedCells; + cells.AsCellSet(castedCells); + auto dims = castedCells.GetSchedulingRange(vtkm::TopologyElementTagPoint()); + vtkm::Vec3f spacing = { static_cast(bounds.X.Length()) / (dims[0] - 1), + static_cast(bounds.Y.Length()) / (dims[1] - 1), + static_cast(bounds.Z.Length()) / (dims[2] - 1) }; + std::cout << spacing << std::endl; + constexpr static vtkm::FloatDefault SPEED_OF_LIGHT = + static_cast(2.99792458e8); + spacing = spacing * spacing; + + + vtkm::Id steps = 50; + vtkm::FloatDefault length = static_cast( + 1.0 / (SPEED_OF_LIGHT * vtkm::Sqrt(1. / spacing[0] + 1. / spacing[1] + 1. / spacing[2]))); + std::cout << "CFL length : " << length << std::endl; + + vtkm::filter::flow::Streamline streamline; + + streamline.SetStepSize(length); + streamline.SetNumberOfSteps(steps); + streamline.SetSeeds(seeds); + streamline.SetVectorFieldType(vtkm::filter::flow::VectorFieldType::ELECTRO_MAGNETIC_FIELD_TYPE); + streamline.SetEField("E"); + streamline.SetBField("B"); + std::cout << "[pre] Executing test" << std::endl; + auto output = streamline.Execute(dataset); + std::cout << "[post] Executing test" << std::endl; + std::cout << "[pre] Executing asserts" << std::endl; + VTKM_TEST_ASSERT(output.GetNumberOfCoordinateSystems() == 1, + "Wrong number of coordinate systems in the output dataset"); + VTKM_TEST_ASSERT(output.GetCoordinateSystem().GetNumberOfPoints() == 2550, + "Wrong number of coordinates"); + VTKM_TEST_ASSERT(output.GetCellSet().GetNumberOfCells() == 50, "Wrong number of cells"); + std::cout << "[post] Executing asserts" << std::endl; +} +} + +int UnitTestStreamlineFilterWarpX(int argc, char* argv[]) +{ + return vtkm::cont::testing::Testing::Run(TestStreamlineFilters, argc, argv); +} -- GitLab From 2004177dbeb825361f02051a7782371563453d4b Mon Sep 17 00:00:00 2001 From: Abhishek Yenpure Date: Mon, 22 Aug 2022 09:21:47 -0700 Subject: [PATCH 3/4] Making test data in control environment --- .../testing/UnitTestStreamlineFilterWarpX.cxx | 54 +++++++------------ 1 file changed, 19 insertions(+), 35 deletions(-) diff --git a/vtkm/filter/flow/testing/UnitTestStreamlineFilterWarpX.cxx b/vtkm/filter/flow/testing/UnitTestStreamlineFilterWarpX.cxx index 7813cb90b..cf4d08119 100644 --- a/vtkm/filter/flow/testing/UnitTestStreamlineFilterWarpX.cxx +++ b/vtkm/filter/flow/testing/UnitTestStreamlineFilterWarpX.cxx @@ -19,33 +19,6 @@ namespace { -class GetChargedParticles : public vtkm::worklet::WorkletMapField -{ -public: - GetChargedParticles() {} - - using ControlSignature = void(FieldIn pos, - FieldIn mom, - FieldIn mass, - FieldIn charge, - FieldIn weighting, - FieldOut electrons); - - using ExecutionSignature = void(WorkIndex, _1, _2, _3, _4, _5, _6); - - void operator()(const vtkm::Id index, - const vtkm::Vec3f& pos, - const vtkm::Vec3f& mom, - const vtkm::Float64& mass, - const vtkm::Float64& charge, - const vtkm::Float64& w, - vtkm::ChargedParticle& electron) const - { - electron = vtkm::ChargedParticle(pos, index, mass, charge, w, mom); - } -}; - - void GenerateChargedParticles(const vtkm::cont::ArrayHandle& pos, const vtkm::cont::ArrayHandle& mom, const vtkm::cont::ArrayHandle& mass, @@ -53,10 +26,24 @@ void GenerateChargedParticles(const vtkm::cont::ArrayHandle& pos, const vtkm::cont::ArrayHandle& weight, vtkm::cont::ArrayHandle& seeds) { - vtkm::cont::Invoker invoker; - invoker(GetChargedParticles{}, pos, mom, mass, charge, weight, seeds); -} + 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, electron); + } +} void TestStreamlineFilters() { @@ -98,7 +85,6 @@ void TestStreamlineFilters() static_cast(2.99792458e8); spacing = spacing * spacing; - vtkm::Id steps = 50; vtkm::FloatDefault length = static_cast( 1.0 / (SPEED_OF_LIGHT * vtkm::Sqrt(1. / spacing[0] + 1. / spacing[1] + 1. / spacing[2]))); @@ -112,16 +98,14 @@ void TestStreamlineFilters() streamline.SetVectorFieldType(vtkm::filter::flow::VectorFieldType::ELECTRO_MAGNETIC_FIELD_TYPE); streamline.SetEField("E"); streamline.SetBField("B"); - std::cout << "[pre] Executing test" << std::endl; + auto output = streamline.Execute(dataset); - std::cout << "[post] Executing test" << std::endl; - std::cout << "[pre] Executing asserts" << std::endl; + VTKM_TEST_ASSERT(output.GetNumberOfCoordinateSystems() == 1, "Wrong number of coordinate systems in the output dataset"); VTKM_TEST_ASSERT(output.GetCoordinateSystem().GetNumberOfPoints() == 2550, "Wrong number of coordinates"); VTKM_TEST_ASSERT(output.GetCellSet().GetNumberOfCells() == 50, "Wrong number of cells"); - std::cout << "[post] Executing asserts" << std::endl; } } -- GitLab From 81b85e8ec1c2adc6e5a29610a98eb76184db3f86 Mon Sep 17 00:00:00 2001 From: Abhishek Yenpure Date: Mon, 22 Aug 2022 10:37:54 -0700 Subject: [PATCH 4/4] Removing unnecessary includes --- vtkm/filter/flow/testing/UnitTestStreamlineFilterWarpX.cxx | 2 -- 1 file changed, 2 deletions(-) diff --git a/vtkm/filter/flow/testing/UnitTestStreamlineFilterWarpX.cxx b/vtkm/filter/flow/testing/UnitTestStreamlineFilterWarpX.cxx index cf4d08119..83b5e66b0 100644 --- a/vtkm/filter/flow/testing/UnitTestStreamlineFilterWarpX.cxx +++ b/vtkm/filter/flow/testing/UnitTestStreamlineFilterWarpX.cxx @@ -10,11 +10,9 @@ #include -#include #include #include #include -#include namespace { -- GitLab