Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-08 13:23:51 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Adding specialization for Particle/Field for XGC

Browse Source 
Adding XGC Field
Adding updates to XGCField
Adding Updates for generalization
Adding WarpXStreamlines and Streamsurface
Adding changes to support XGC Poincare
Finalizing XGC analysis
This commit is contained in:
Abhishek Yenpure 2022-09-25 12:42:23 -07:00
parent 34f07c37cf
commit e6291a9b91
46 changed files with 3202 additions and 1270 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

37
vtkm/Particle.h
View File

@ -97,11 +97,13 @@ public:
Particle(const vtkm::Vec3f& p,
const vtkm::Id& id,
const vtkm::Id& numSteps = 0,
const vtkm::Id& numPunctures = 0,
const vtkm::ParticleStatus& status = vtkm::ParticleStatus(),
const vtkm::FloatDefault& time = 0)
: Position(p)
, ID(id)
, NumSteps(numSteps)
, NumPunctures(numPunctures)
, Status(status)
, Time(time)
{
@ -112,6 +114,7 @@ public:
: Position(p.Position)
, ID(p.ID)
, NumSteps(p.NumSteps)
, NumPunctures(p.NumPunctures)
, Status(p.Status)
, Time(p.Time)
{
@ -134,6 +137,9 @@ public:
VTKM_EXEC_CONT vtkm::Id GetNumberOfSteps() const { return this->NumSteps; }
VTKM_EXEC_CONT void SetNumberOfSteps(vtkm::Id numSteps) { this->NumSteps = numSteps; }
VTKM_EXEC_CONT vtkm::Id GetNumberOfPunctures() const { return this->NumPunctures; }
VTKM_EXEC_CONT void SetNumberOfPunctures(vtkm::Id punctures) { this->NumPunctures = punctures; }
VTKM_EXEC_CONT vtkm::ParticleStatus GetStatus() const { return this->Status; }
VTKM_EXEC_CONT vtkm::ParticleStatus& GetStatus() { return this->Status; }
VTKM_EXEC_CONT void SetStatus(vtkm::ParticleStatus status) { this->Status = status; }
@ -143,7 +149,7 @@ public:
VTKM_EXEC_CONT
vtkm::Vec3f Velocity(const vtkm::VecVariable<vtkm::Vec3f, 2>& vectors,
const vtkm::FloatDefault& vtkmNotUsed(length))
const vtkm::FloatDefault& vtkmNotUsed(length)) const
{
// Velocity is evaluated from the Velocity field
// and is not influenced by the particle
@ -169,6 +175,7 @@ private:
vtkm::Vec3f Position;
vtkm::Id ID = -1;
vtkm::Id NumSteps = 0;
vtkm::Id NumPunctures = 0;
vtkm::ParticleStatus Status;
vtkm::FloatDefault Time = 0;
@ -178,6 +185,7 @@ public:
constexpr std::size_t sz = sizeof(vtkm::Vec3f) // Pos
+ sizeof(vtkm::Id) // ID
+ sizeof(vtkm::Id) // NumSteps
+ sizeof(vtkm::Id) // NumPunctures
+ sizeof(vtkm::UInt8) // Status
+ sizeof(vtkm::FloatDefault); // Time
@ -244,7 +252,7 @@ public:
VTKM_EXEC_CONT
vtkm::Vec3f Velocity(const vtkm::VecVariable<vtkm::Vec3f, 2>& vectors,
const vtkm::FloatDefault& length)
const vtkm::FloatDefault& length) const
{
VTKM_ASSERT(vectors.GetNumberOfComponents() == 2);
@ -302,7 +310,7 @@ private:
vtkm::Float64 Mass;
vtkm::Float64 Charge;
vtkm::Float64 Weighting;
vtkm::Vec3f Momentum;
mutable vtkm::Vec3f Momentum;
constexpr static vtkm::FloatDefault SPEED_OF_LIGHT =
static_cast<vtkm::FloatDefault>(2.99792458e8);
@ -366,7 +374,30 @@ public:
p.SetTime(time);
}
};
/*
template <>
struct Serialization<vtkm::XGCParticle>
{
public:
static VTKM_CONT void save(BinaryBuffer& bb, const vtkm::XGCParticle& p)
{
vtkmdiy::save(bb, p.Pos);
vtkmdiy::save(bb, p.ID);
vtkmdiy::save(bb, p.NumSteps);
vtkmdiy::save(bb, p.Status);
vtkmdiy::save(bb, p.Time);
}
static VTKM_CONT void load(BinaryBuffer& bb, vtkm::XGCParticle& p)
{
vtkmdiy::load(bb, p.Pos);
vtkmdiy::load(bb, p.ID);
vtkmdiy::load(bb, p.NumSteps);
vtkmdiy::load(bb, p.Status);
vtkmdiy::load(bb, p.Time);
}
};
*/
template <>
struct Serialization<vtkm::ChargedParticle>
{

14
vtkm/filter/flow/CMakeLists.txt
View File

@ -20,15 +20,13 @@ set(flow_headers
PathParticle.h
Streamline.h
StreamSurface.h
WarpXStreamline.h
XGCPoincare.h
)
set(flow_sources
internal/Messenger.cxx
FilterParticleAdvection.cxx
ParticleAdvection.cxx
Pathline.cxx
PathParticle.cxx
Streamline.cxx
)
set(flow_device_sources
@ -37,6 +35,14 @@ set(flow_device_sources
FilterParticleAdvectionSteadyState.cxx
FilterParticleAdvectionUnsteadyState.cxx
StreamSurface.cxx
Lagrangian.cxx
LagrangianStructures.cxx
ParticleAdvection.cxx
Pathline.cxx
PathParticle.cxx
Streamline.cxx
WarpXStreamline.cxx
XGCPoincare.cxx
)
vtkm_library(

19
vtkm/filter/flow/FilterParticleAdvection.h
View File

@ -68,23 +68,6 @@ public:
this->SolverType = vtkm::filter::flow::IntegrationSolverType::EULER_TYPE;
}
VTKM_CONT
void SetVectorFieldType(vtkm::filter::flow::VectorFieldType vecFieldType)
{
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; }
@ -104,10 +87,10 @@ public:
protected:
VTKM_CONT virtual void ValidateOptions() const;
VTKM_CONT virtual vtkm::filter::flow::FlowResultType GetResultType() const = 0;
bool BlockIdsSet = false;
std::vector<vtkm::Id> BlockIds;
vtkm::Id NumberOfSteps = 0;
vtkm::cont::UnknownArrayHandle Seeds;
vtkm::filter::flow::IntegrationSolverType SolverType =

85
vtkm/filter/flow/FilterParticleAdvectionSteadyState.cxx
View File

@ -9,10 +9,6 @@
//============================================================================
#include <vtkm/filter/flow/FilterParticleAdvectionSteadyState.h>
#include <vtkm/filter/flow/internal/DataSetIntegratorSteadyState.h>
#include <vtkm/filter/flow/internal/BoundsMap.h>
#include <vtkm/filter/flow/internal/ParticleAdvector.h>
namespace vtkm
{
@ -21,87 +17,6 @@ namespace filter
namespace flow
{
namespace
{
VTKM_CONT std::vector<vtkm::filter::flow::internal::DataSetIntegratorSteadyState>
CreateDataSetIntegrators(
const vtkm::cont::PartitionedDataSet& input,
const vtkm::cont::Variant<std::string, std::pair<std::string, std::string>>& 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;
std::vector<DSIType> dsi;
for (vtkm::Id i = 0; i < input.GetNumberOfPartitions(); i++)
{
vtkm::Id blockId = boundsMap.GetLocalBlockId(i);
auto ds = input.GetPartition(i);
if (activeField.IsType<DSIType::VelocityFieldNameType>())
{
const auto& fieldNm = activeField.Get<DSIType::VelocityFieldNameType>();
if (!ds.HasPointField(fieldNm) && !ds.HasCellField(fieldNm))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
}
else if (activeField.IsType<DSIType::ElectroMagneticFieldNameType>())
{
const auto& fieldNm = activeField.Get<DSIType::ElectroMagneticFieldNameType>();
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;
}
} //anonymous namespace
VTKM_CONT vtkm::cont::PartitionedDataSet FilterParticleAdvectionSteadyState::DoExecutePartitions(
const vtkm::cont::PartitionedDataSet& input)
{
using DSIType = vtkm::filter::flow::internal::DataSetIntegratorSteadyState;
this->ValidateOptions();
using VariantType = vtkm::cont::Variant<std::string, std::pair<std::string, std::string>>;
VariantType variant;
if (this->VecFieldType == vtkm::filter::flow::VectorFieldType::VELOCITY_FIELD_TYPE)
{
const auto& field = this->GetActiveFieldName();
variant.Emplace<DSIType::VelocityFieldNameType>(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<DSIType::ElectroMagneticFieldNameType>(electric, magnetic);
}
vtkm::filter::flow::internal::BoundsMap boundsMap;
if (this->BlockIdsSet)
boundsMap = vtkm::filter::flow::internal::BoundsMap(input, this->BlockIds);
else
boundsMap = vtkm::filter::flow::internal::BoundsMap(input);
auto dsi = CreateDataSetIntegrators(
input, variant, boundsMap, this->SolverType, this->VecFieldType, this->GetResultType());
vtkm::filter::flow::internal::ParticleAdvector<DSIType> pav(boundsMap,
dsi,
this->UseThreadedAlgorithm,
this->UseAsynchronousCommunication,
this->GetResultType());
return pav.Execute(this->NumberOfSteps, this->StepSize, this->Seeds);
}
}
}
} // namespace vtkm::filter::flow

70
vtkm/filter/flow/FilterParticleAdvectionSteadyState.h
View File

@ -15,19 +15,81 @@
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
#include <vtkm/filter/flow/internal/BoundsMap.h>
#include <vtkm/filter/flow/internal/DataSetIntegratorSteadyState.h>
#include <vtkm/filter/flow/internal/ParticleAdvector.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
template <typename Derived>
struct FlowTraits;
template <typename Derived>
class VTKM_FILTER_FLOW_EXPORT FilterParticleAdvectionSteadyState : public FilterParticleAdvection
{
private:
VTKM_CONT vtkm::cont::PartitionedDataSet DoExecutePartitions(
const vtkm::cont::PartitionedDataSet& inData) override;
};
public:
using ParticleType = typename FlowTraits<Derived>::ParticleType;
using FieldType = typename FlowTraits<Derived>::FieldType;
using TerminationType = typename FlowTraits<Derived>::TerminationType;
using AnalysisType = typename FlowTraits<Derived>::AnalysisType;
private:
VTKM_CONT FieldType GetField(const vtkm::cont::DataSet& data) const
{
const Derived* inst = static_cast<const Derived*>(this);
return inst->GetField(data);
}
VTKM_CONT TerminationType GetTermination(const vtkm::cont::DataSet& data) const
{
const Derived* inst = static_cast<const Derived*>(this);
return inst->GetTermination(data);
}
VTKM_CONT AnalysisType GetAnalysis(const vtkm::cont::DataSet& data) const
{
const Derived* inst = static_cast<const Derived*>(this);
return inst->GetAnalysis(data);
}
VTKM_CONT vtkm::cont::PartitionedDataSet DoExecutePartitions(
const vtkm::cont::PartitionedDataSet& input)
{
this->ValidateOptions();
using DSIType = vtkm::filter::flow::internal::
DataSetIntegratorSteadyState<ParticleType, FieldType, TerminationType, AnalysisType>;
vtkm::filter::flow::internal::BoundsMap boundsMap(input);
std::vector<DSIType> dsi;
for (vtkm::Id i = 0; i < input.GetNumberOfPartitions(); i++)
{
vtkm::Id blockId = boundsMap.GetLocalBlockId(i);
auto dataset = input.GetPartition(i);
// Build the field for the current dataset
FieldType field = this->GetField(dataset);
// Build the termination for the current dataset
TerminationType termination = this->GetTermination(dataset);
// Build the analysis for the current dataset
AnalysisType analysis = this->GetAnalysis(dataset);
dsi.emplace_back(blockId, field, dataset, this->SolverType, termination, analysis);
}
vtkm::filter::flow::internal::ParticleAdvector<DSIType> pav(
boundsMap, dsi, this->UseThreadedAlgorithm);
vtkm::cont::ArrayHandle<ParticleType> particles;
this->Seeds.AsArrayHandle(particles);
return pav.Execute(particles, this->StepSize);
}
};
}
}
} // namespace vtkm::filter::flow

67
vtkm/filter/flow/FilterParticleAdvectionUnsteadyState.cxx
View File

@ -9,8 +9,6 @@
//============================================================================
#include <vtkm/filter/flow/FilterParticleAdvectionUnsteadyState.h>
#include <vtkm/filter/flow/internal/DataSetIntegratorUnsteadyState.h>
#include <vtkm/filter/flow/internal/ParticleAdvector.h>
namespace vtkm
{
@ -18,71 +16,6 @@ namespace filter
{
namespace flow
{
namespace
{
VTKM_CONT std::vector<vtkm::filter::flow::internal::DataSetIntegratorUnsteadyState>
CreateDataSetIntegrators(const vtkm::cont::PartitionedDataSet& input,
const vtkm::cont::PartitionedDataSet& input2,
const std::string& activeField,
const vtkm::FloatDefault timer1,
const vtkm::FloatDefault timer2,
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::DataSetIntegratorUnsteadyState;
std::vector<DSIType> dsi;
for (vtkm::Id i = 0; i < input.GetNumberOfPartitions(); i++)
{
vtkm::Id blockId = boundsMap.GetLocalBlockId(i);
auto ds1 = input.GetPartition(i);
auto ds2 = input2.GetPartition(i);
if ((!ds1.HasPointField(activeField) && !ds1.HasCellField(activeField)) ||
(!ds2.HasPointField(activeField) && !ds2.HasCellField(activeField)))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
dsi.emplace_back(
ds1, ds2, timer1, timer2, blockId, activeField, solverType, vecFieldType, resultType);
}
return dsi;
}
} // anonymous namespace
VTKM_CONT void FilterParticleAdvectionUnsteadyState::ValidateOptions() const
{
this->FilterParticleAdvection::ValidateOptions();
if (this->Time1 >= this->Time2)
throw vtkm::cont::ErrorFilterExecution("PreviousTime must be less than NextTime");
}
VTKM_CONT vtkm::cont::PartitionedDataSet FilterParticleAdvectionUnsteadyState::DoExecutePartitions(
const vtkm::cont::PartitionedDataSet& input)
{
using DSIType = vtkm::filter::flow::internal::DataSetIntegratorUnsteadyState;
this->ValidateOptions();
vtkm::filter::flow::internal::BoundsMap boundsMap(input);
auto dsi = CreateDataSetIntegrators(input,
this->Input2,
this->GetActiveFieldName(),
this->Time1,
this->Time2,
boundsMap,
this->SolverType,
this->VecFieldType,
this->GetResultType());
vtkm::filter::flow::internal::ParticleAdvector<DSIType> pav(boundsMap,
dsi,
this->UseThreadedAlgorithm,
this->UseAsynchronousCommunication,
this->GetResultType());
return pav.Execute(this->NumberOfSteps, this->StepSize, this->Seeds);
}
}
}

78
vtkm/filter/flow/FilterParticleAdvectionUnsteadyState.h
View File

@ -14,6 +14,10 @@
#include <vtkm/filter/flow/FilterParticleAdvection.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
#include <vtkm/filter/flow/internal/BoundsMap.h>
#include <vtkm/filter/flow/internal/DataSetIntegratorUnsteadyState.h>
#include <vtkm/filter/flow/internal/ParticleAdvector.h>
namespace vtkm
{
namespace filter
@ -21,9 +25,18 @@ namespace filter
namespace flow
{
template <typename Derived>
struct FlowTraits;
template <typename Derived>
class VTKM_FILTER_FLOW_EXPORT FilterParticleAdvectionUnsteadyState : public FilterParticleAdvection
{
public:
using ParticleType = typename FlowTraits<Derived>::ParticleType;
using FieldType = typename FlowTraits<Derived>::FieldType;
using TerminationType = typename FlowTraits<Derived>::TerminationType;
using AnalysisType = typename FlowTraits<Derived>::AnalysisType;
VTKM_CONT void SetPreviousTime(vtkm::FloatDefault t1) { this->Time1 = t1; }
VTKM_CONT void SetNextTime(vtkm::FloatDefault t2) { this->Time2 = t2; }
@ -32,12 +45,69 @@ public:
VTKM_CONT void SetNextDataSet(const vtkm::cont::PartitionedDataSet& pds) { this->Input2 = pds; }
protected:
VTKM_CONT virtual void ValidateOptions() const override;
private:
VTKM_CONT FieldType GetField(const vtkm::cont::DataSet& data) const
{
const Derived* inst = static_cast<const Derived*>(this);
return inst->GetField(data);
}
VTKM_CONT TerminationType GetTermination(const vtkm::cont::DataSet& data) const
{
const Derived* inst = static_cast<const Derived*>(this);
return inst->GetTermination(data);
}
VTKM_CONT AnalysisType GetAnalysis(const vtkm::cont::DataSet& data) const
{
const Derived* inst = static_cast<const Derived*>(this);
return inst->GetAnalysis(data);
}
VTKM_CONT vtkm::cont::PartitionedDataSet DoExecutePartitions(
const vtkm::cont::PartitionedDataSet& inData) override;
const vtkm::cont::PartitionedDataSet& input)
{
this->ValidateOptions();
using DSIType = vtkm::filter::flow::internal::
DataSetIntegratorUnsteadyState<ParticleType, FieldType, TerminationType, AnalysisType>;
vtkm::filter::flow::internal::BoundsMap boundsMap(input);
std::vector<DSIType> dsi;
for (vtkm::Id i = 0; i < input.GetNumberOfPartitions(); i++)
{
vtkm::Id blockId = boundsMap.GetLocalBlockId(i);
auto ds1 = input.GetPartition(i);
auto ds2 = this->Input2.GetPartition(i);
// Build the field for the current dataset
FieldType field1 = this->GetField(ds1);
FieldType field2 = this->GetField(ds2);
// Build the termination for the current dataset
TerminationType termination = this->GetTermination(ds1);
AnalysisType analysis = this->GetAnalysis(ds1);
dsi.emplace_back(blockId,
field1,
field2,
ds1,
ds2,
this->Time1,
this->Time2,
this->SolverType,
termination,
analysis);
}
vtkm::filter::flow::internal::ParticleAdvector<DSIType> pav(
boundsMap, dsi, this->UseThreadedAlgorithm);
vtkm::cont::ArrayHandle<ParticleType> particles;
this->Seeds.AsArrayHandle(particles);
return pav.Execute(particles, this->StepSize);
}
vtkm::cont::PartitionedDataSet Input2;
vtkm::FloatDefault Time1 = -1;

16
vtkm/filter/flow/Lagrangian.cxx
View File

@ -21,6 +21,14 @@
#include <vtkm/filter/flow/worklet/RK4Integrator.h>
#include <vtkm/filter/flow/worklet/Stepper.h>
#include <vtkm/filter/flow/worklet/Analysis.h>
#include <vtkm/filter/flow/worklet/Termination.h>
//#include <vtkm/worklet/WorkletMapField.h>
//#include <cstring>
//#include <sstream>
//#include <string.h>
namespace vtkm
{
namespace filter
@ -220,7 +228,6 @@ VTKM_CONT vtkm::cont::DataSet Lagrangian::DoExecute(const vtkm::cont::DataSet& i
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
vtkm::worklet::flow::ParticleAdvection particleadvection;
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> res;
const auto field = input.GetField(this->GetActiveFieldName());
FieldType velocities(field.GetData().AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Vec3f>>(),
@ -228,9 +235,10 @@ VTKM_CONT vtkm::cont::DataSet Lagrangian::DoExecute(const vtkm::cont::DataSet& i
GridEvalType gridEval(coords, cells, velocities);
Stepper rk4(gridEval, static_cast<vtkm::Float32>(this->StepSize));
res = particleadvection.Run(rk4, basisParticleArray, 1); // Taking a single step
auto particles = res.Particles;
vtkm::worklet::flow::NormalTermination termination(1);
vtkm::worklet::flow::NoAnalysis<vtkm::Particle> analysis;
particleadvection.Run(rk4, basisParticleArray, termination, analysis); // Taking a single step
auto particles = analysis.Particles;
vtkm::cont::DataSet outputData;

14
vtkm/filter/flow/LagrangianStructures.cxx
View File

@ -13,14 +13,17 @@
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/DataSetBuilderUniform.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/cont/Invoker.h>
#include <vtkm/filter/flow/LagrangianStructures.h>
#include <vtkm/filter/flow/worklet/Analysis.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/GridEvaluators.h>
#include <vtkm/filter/flow/worklet/LagrangianStructures.h>
#include <vtkm/filter/flow/worklet/ParticleAdvection.h>
#include <vtkm/filter/flow/worklet/RK4Integrator.h>
#include <vtkm/filter/flow/worklet/Stepper.h>
#include <vtkm/filter/flow/worklet/Termination.h>
namespace
{
@ -144,11 +147,14 @@ VTKM_CONT vtkm::cont::DataSet LagrangianStructures::DoExecute(const vtkm::cont::
GridEvaluator evaluator(input.GetCoordinateSystem(), input.GetCellSet(), velocities);
Stepper integrator(evaluator, stepSize);
vtkm::worklet::flow::ParticleAdvection particles;
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> advectionResult;
vtkm::worklet::flow::NormalTermination termination(numberOfSteps);
vtkm::worklet::flow::NoAnalysis<vtkm::Particle> analysis;
vtkm::cont::ArrayHandle<vtkm::Particle> advectionPoints;
this->Invoke(detail::MakeParticles{}, lcsInputPoints, advectionPoints);
advectionResult = particles.Run(integrator, advectionPoints, numberOfSteps);
this->Invoke(detail::ExtractParticlePosition{}, advectionResult.Particles, lcsOutputPoints);
vtkm::cont::Invoker invoke;
invoke(detail::MakeParticles{}, lcsInputPoints, advectionPoints);
particles.Run(integrator, advectionPoints, termination, analysis);
invoke(detail::ExtractParticlePosition{}, analysis.Particles, lcsOutputPoints);
}
// FTLE output field
vtkm::cont::ArrayHandle<vtkm::FloatDefault> outputField;

33
vtkm/filter/flow/ParticleAdvection.cxx
View File

@ -17,9 +17,38 @@ namespace filter
namespace flow
{
VTKM_CONT vtkm::filter::flow::FlowResultType ParticleAdvection::GetResultType() const
// using ParticleType = vtkm::Particle;
// using TerminationType = vtkm::worklet::flow::NormalTermination;
// using AnalysisType = vtkm::worklet::flow::Streamline;
// using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
// using FieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
VTKM_CONT ParticleAdvection::FieldType ParticleAdvection::GetField(
const vtkm::cont::DataSet& dataset) const
{
return vtkm::filter::flow::FlowResultType::PARTICLE_ADVECT_TYPE;
const auto& fieldNm = this->GetActiveFieldName();
if (!dataset.HasPointField(fieldNm) && !dataset.HasCellField(fieldNm))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
auto assoc = dataset.GetField(fieldNm).GetAssociation();
ArrayType arr;
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(fieldNm).GetData(), arr);
return ParticleAdvection::FieldType(arr, assoc);
}
VTKM_CONT ParticleAdvection::TerminationType ParticleAdvection::GetTermination(
const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return ParticleAdvection::TerminationType(this->NumberOfSteps);
}
VTKM_CONT ParticleAdvection::AnalysisType ParticleAdvection::GetAnalysis(
const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return ParticleAdvection::AnalysisType{};
}
}

31
vtkm/filter/flow/ParticleAdvection.h
View File

@ -15,6 +15,10 @@
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
#include <vtkm/filter/flow/worklet/Analysis.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/Termination.h>
namespace vtkm
{
namespace filter
@ -22,17 +26,38 @@ namespace filter
namespace flow
{
class ParticleAdvection;
template <>
struct FlowTraits<ParticleAdvection>
{
using ParticleType = vtkm::Particle;
using TerminationType = vtkm::worklet::flow::NormalTermination;
using AnalysisType = vtkm::worklet::flow::NoAnalysis<ParticleType>;
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
};
/// \brief Advect particles in a vector field.
/// Takes as input a vector field and seed locations and generates the
/// end points for each seed through the vector field.
class VTKM_FILTER_FLOW_EXPORT ParticleAdvection
: public vtkm::filter::flow::FilterParticleAdvectionSteadyState
: public vtkm::filter::flow::FilterParticleAdvectionSteadyState<ParticleAdvection>
{
public:
private:
VTKM_CONT vtkm::filter::flow::FlowResultType GetResultType() const override;
using ParticleType = typename FlowTraits<ParticleAdvection>::ParticleType;
using TerminationType = typename FlowTraits<ParticleAdvection>::TerminationType;
using AnalysisType = typename FlowTraits<ParticleAdvection>::AnalysisType;
using ArrayType = typename FlowTraits<ParticleAdvection>::ArrayType;
using FieldType = typename FlowTraits<ParticleAdvection>::FieldType;
VTKM_CONT FieldType GetField(const vtkm::cont::DataSet& data) const;
VTKM_CONT TerminationType GetTermination(const vtkm::cont::DataSet& data) const;
VTKM_CONT AnalysisType GetAnalysis(const vtkm::cont::DataSet& data) const;
};
}

31
vtkm/filter/flow/PathParticle.cxx
View File

@ -17,11 +17,38 @@ namespace filter
namespace flow
{
VTKM_CONT vtkm::filter::flow::FlowResultType PathParticle::GetResultType() const
VTKM_CONT PathParticle::FieldType PathParticle::GetField(const vtkm::cont::DataSet& dataset) const
{
return vtkm::filter::flow::FlowResultType::PARTICLE_ADVECT_TYPE;
const auto& fieldNm = this->GetActiveFieldName();
if (!dataset.HasPointField(fieldNm) && !dataset.HasCellField(fieldNm))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
auto assoc = dataset.GetField(fieldNm).GetAssociation();
ArrayType arr;
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(fieldNm).GetData(), arr);
return PathParticle::FieldType(arr, assoc);
}
VTKM_CONT PathParticle::TerminationType PathParticle::GetTermination(
const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return PathParticle::TerminationType(this->NumberOfSteps);
}
VTKM_CONT PathParticle::AnalysisType PathParticle::GetAnalysis(
const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return PathParticle::AnalysisType();
}
//VTKM_CONT vtkm::filter::flow::FlowResultType PathParticle::GetResultType() const
//{
// return vtkm::filter::flow::FlowResultType::PARTICLE_ADVECT_TYPE;
//}
}
}
} // namespace vtkm::filter::flow

32
vtkm/filter/flow/PathParticle.h
View File

@ -15,6 +15,10 @@
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
#include <vtkm/filter/flow/worklet/Analysis.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/Termination.h>
namespace vtkm
{
namespace filter
@ -22,16 +26,38 @@ namespace filter
namespace flow
{
class PathParticle;
template <>
struct FlowTraits<PathParticle>
{
using ParticleType = vtkm::Particle;
using TerminationType = vtkm::worklet::flow::NormalTermination;
using AnalysisType = vtkm::worklet::flow::NoAnalysis<ParticleType>;
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
};
/// \brief Advect particles in a vector field.
/// Takes as input a vector field and seed locations and generates the
/// end points for each seed through the vector field.
class VTKM_FILTER_FLOW_EXPORT PathParticle
: public vtkm::filter::flow::FilterParticleAdvectionUnsteadyState
: public vtkm::filter::flow::FilterParticleAdvectionUnsteadyState<PathParticle>
{
private:
VTKM_CONT vtkm::filter::flow::FlowResultType GetResultType() const override;
public:
using ParticleType = typename FlowTraits<PathParticle>::ParticleType;
using TerminationType = typename FlowTraits<PathParticle>::TerminationType;
using AnalysisType = typename FlowTraits<PathParticle>::AnalysisType;
using ArrayType = typename FlowTraits<PathParticle>::ArrayType;
using FieldType = typename FlowTraits<PathParticle>::FieldType;
VTKM_CONT FieldType GetField(const vtkm::cont::DataSet& data) const;
VTKM_CONT TerminationType GetTermination(const vtkm::cont::DataSet& data) const;
VTKM_CONT AnalysisType GetAnalysis(const vtkm::cont::DataSet& data) const;
};
}

29
vtkm/filter/flow/Pathline.cxx
View File

@ -17,11 +17,36 @@ namespace filter
namespace flow
{
VTKM_CONT vtkm::filter::flow::FlowResultType Pathline::GetResultType() const
VTKM_CONT Pathline::FieldType Pathline::GetField(const vtkm::cont::DataSet& dataset) const
{
return vtkm::filter::flow::FlowResultType::STREAMLINE_TYPE;
const auto& fieldNm = this->GetActiveFieldName();
if (!dataset.HasPointField(fieldNm) && !dataset.HasCellField(fieldNm))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
auto assoc = dataset.GetField(fieldNm).GetAssociation();
ArrayType arr;
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(fieldNm).GetData(), arr);
return vtkm::worklet::flow::VelocityField<ArrayType>(arr, assoc);
}
VTKM_CONT Pathline::TerminationType Pathline::GetTermination(
const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return Pathline::TerminationType(this->NumberOfSteps);
}
VTKM_CONT Pathline::AnalysisType Pathline::GetAnalysis(const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return Pathline::AnalysisType(this->NumberOfSteps);
}
//VTKM_CONT vtkm::filter::flow::FlowResultType Pathline::GetResultType() const
//{
// return vtkm::filter::flow::FlowResultType::STREAMLINE_TYPE;
//}
}
}
} // namespace vtkm::filter::flow

32
vtkm/filter/flow/Pathline.h
View File

@ -15,6 +15,10 @@
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
#include <vtkm/filter/flow/worklet/Analysis.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/Termination.h>
namespace vtkm
{
namespace filter
@ -22,16 +26,38 @@ namespace filter
namespace flow
{
class Pathline;
template <>
struct FlowTraits<Pathline>
{
using ParticleType = vtkm::Particle;
using TerminationType = vtkm::worklet::flow::NormalTermination;
using AnalysisType = vtkm::worklet::flow::StreamlineAnalysis<ParticleType>;
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
};
/// \brief Advect particles in a vector field.
/// Takes as input a vector field and seed locations and generates the
/// end points for each seed through the vector field.
class VTKM_FILTER_FLOW_EXPORT Pathline
: public vtkm::filter::flow::FilterParticleAdvectionUnsteadyState
: public vtkm::filter::flow::FilterParticleAdvectionUnsteadyState<Pathline>
{
private:
VTKM_CONT vtkm::filter::flow::FlowResultType GetResultType() const override;
public:
using ParticleType = typename FlowTraits<Pathline>::ParticleType;
using TerminationType = typename FlowTraits<Pathline>::TerminationType;
using AnalysisType = typename FlowTraits<Pathline>::AnalysisType;
using ArrayType = typename FlowTraits<Pathline>::ArrayType;
using FieldType = typename FlowTraits<Pathline>::FieldType;
VTKM_CONT FieldType GetField(const vtkm::cont::DataSet& data) const;
VTKM_CONT TerminationType GetTermination(const vtkm::cont::DataSet& data) const;
VTKM_CONT AnalysisType GetAnalysis(const vtkm::cont::DataSet& data) const;
};
}

19
vtkm/filter/flow/StreamSurface.cxx
View File

@ -1,4 +1,5 @@
//============================================================================
//=============================================================================
//
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
@ -6,7 +7,8 @@
// 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 <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ErrorFilterExecution.h>
@ -65,19 +67,22 @@ VTKM_CONT vtkm::cont::DataSet StreamSurface::DoExecute(const vtkm::cont::DataSet
GridEvalType eval(coords, cells, velocities);
Stepper rk4(eval, this->StepSize);
vtkm::worklet::flow::Streamline streamline;
using ParticleArray = vtkm::cont::ArrayHandle<vtkm::Particle>;
vtkm::cont::ArrayHandle<vtkm::Particle> seedArray;
vtkm::cont::ArrayCopy(this->Seeds.AsArrayHandle<ParticleArray>(), seedArray);
auto res = streamline.Run(rk4, seedArray, this->NumberOfSteps);
vtkm::worklet::flow::ParticleAdvection worklet;
vtkm::worklet::flow::NormalTermination termination(this->NumberOfSteps);
vtkm::worklet::flow::StreamlineAnalysis<vtkm::Particle> analysis(this->NumberOfSteps);
worklet.Run(rk4, seedArray, termination, analysis);
//compute surface from streamlines
vtkm::worklet::flow::StreamSurface streamSurface;
vtkm::cont::ArrayHandle<vtkm::Vec3f> srfPoints;
vtkm::cont::CellSetSingleType<> srfCells;
vtkm::cont::CoordinateSystem slCoords("coordinates", res.Positions);
streamSurface.Run(slCoords, res.PolyLines, srfPoints, srfCells);
vtkm::cont::CoordinateSystem slCoords("coordinates", analysis.Streams);
streamSurface.Run(slCoords, analysis.PolyLines, srfPoints, srfCells);
vtkm::cont::DataSet outData;
vtkm::cont::CoordinateSystem outputCoords("coordinates", srfPoints);

31
vtkm/filter/flow/Streamline.cxx
View File

@ -17,9 +17,36 @@ namespace filter
namespace flow
{
VTKM_CONT vtkm::filter::flow::FlowResultType Streamline::GetResultType() const
// using ParticleType = vtkm::Particle;
// using TerminationType = vtkm::worklet::flow::NormalTermination;
// using AnalysisType = vtkm::worklet::flow::Streamline;
// using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
// using FieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
VTKM_CONT Streamline::FieldType Streamline::GetField(const vtkm::cont::DataSet& dataset) const
{
return vtkm::filter::flow::FlowResultType::STREAMLINE_TYPE;
const auto& fieldNm = this->GetActiveFieldName();
if (!dataset.HasPointField(fieldNm) && !dataset.HasCellField(fieldNm))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
auto assoc = dataset.GetField(fieldNm).GetAssociation();
ArrayType arr;
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(fieldNm).GetData(), arr);
return vtkm::worklet::flow::VelocityField<ArrayType>(arr, assoc);
}
VTKM_CONT Streamline::TerminationType Streamline::GetTermination(
const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return Streamline::TerminationType(this->NumberOfSteps);
}
VTKM_CONT Streamline::AnalysisType Streamline::GetAnalysis(const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return Streamline::AnalysisType(this->NumberOfSteps);
}
}

32
vtkm/filter/flow/Streamline.h
View File

@ -15,6 +15,10 @@
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
#include <vtkm/filter/flow/worklet/Analysis.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/Termination.h>
namespace vtkm
{
namespace filter
@ -22,16 +26,38 @@ namespace filter
namespace flow
{
class Streamline;
template <>
struct FlowTraits<Streamline>
{
using ParticleType = vtkm::Particle;
using TerminationType = vtkm::worklet::flow::NormalTermination;
using AnalysisType = vtkm::worklet::flow::StreamlineAnalysis<ParticleType>;
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
};
/// \brief Advect particles in a vector field.
/// Takes as input a vector field and seed locations and generates the
/// end points for each seed through the vector field.
class VTKM_FILTER_FLOW_EXPORT Streamline
: public vtkm::filter::flow::FilterParticleAdvectionSteadyState
: public vtkm::filter::flow::FilterParticleAdvectionSteadyState<Streamline>
{
private:
VTKM_CONT vtkm::filter::flow::FlowResultType GetResultType() const override;
public:
using ParticleType = typename FlowTraits<Streamline>::ParticleType;
using TerminationType = typename FlowTraits<Streamline>::TerminationType;
using AnalysisType = typename FlowTraits<Streamline>::AnalysisType;
using ArrayType = typename FlowTraits<Streamline>::ArrayType;
using FieldType = typename FlowTraits<Streamline>::FieldType;
VTKM_CONT FieldType GetField(const vtkm::cont::DataSet& data) const;
VTKM_CONT TerminationType GetTermination(const vtkm::cont::DataSet& data) const;
VTKM_CONT AnalysisType GetAnalysis(const vtkm::cont::DataSet& data) const;
};
}

59
vtkm/filter/flow/WarpXStreamline.cxx Normal file
View File

@ -0,0 +1,59 @@
//============================================================================
// 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 <vtkm/filter/flow/WarpXStreamline.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
VTKM_CONT WarpXStreamline::FieldType WarpXStreamline::GetField(
const vtkm::cont::DataSet& dataset) const
{
const auto& electric = this->GetEField();
const auto& magnetic = this->GetBField();
if (!dataset.HasPointField(electric) && !dataset.HasCellField(electric))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
if (!dataset.HasPointField(magnetic) && !dataset.HasCellField(magnetic))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
auto eAssoc = dataset.GetField(electric).GetAssociation();
auto bAssoc = 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(dataset.GetField(electric).GetData(), eField);
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(magnetic).GetData(), bField);
return vtkm::worklet::flow::ElectroMagneticField<ArrayType>(eField, bField, eAssoc);
}
VTKM_CONT WarpXStreamline::TerminationType WarpXStreamline::GetTermination(
const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return WarpXStreamline::TerminationType(this->NumberOfSteps);
}
VTKM_CONT WarpXStreamline::AnalysisType WarpXStreamline::GetAnalysis(
const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return WarpXStreamline::AnalysisType(this->NumberOfSteps);
}
}
}
} // namespace vtkm::filter::flow

74
vtkm/filter/flow/WarpXStreamline.h Normal file
View File

@ -0,0 +1,74 @@
//============================================================================
// 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_filter_flow_WarpXStreamline_h
#define vtk_m_filter_flow_WarpXStreamline_h
#include <vtkm/filter/flow/FilterParticleAdvectionSteadyState.h>
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
#include <vtkm/filter/flow/worklet/Analysis.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/Termination.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
class WarpXStreamline;
template <>
struct FlowTraits<WarpXStreamline>
{
using ParticleType = vtkm::ChargedParticle;
using TerminationType = vtkm::worklet::flow::NormalTermination;
using AnalysisType = vtkm::worklet::flow::StreamlineAnalysis<ParticleType>;
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::flow::ElectroMagneticField<ArrayType>;
};
/// \brief Advect particles in a vector field.
/// Takes as input a vector field and seed locations and generates the
/// end points for each seed through the vector field.
class VTKM_FILTER_FLOW_EXPORT WarpXStreamline
: public vtkm::filter::flow::FilterParticleAdvectionSteadyState<WarpXStreamline>
{
public:
using ParticleType = typename FlowTraits<WarpXStreamline>::ParticleType;
using TerminationType = typename FlowTraits<WarpXStreamline>::TerminationType;
using AnalysisType = typename FlowTraits<WarpXStreamline>::AnalysisType;
using ArrayType = typename FlowTraits<WarpXStreamline>::ArrayType;
using FieldType = typename FlowTraits<WarpXStreamline>::FieldType;
VTKM_CONT FieldType GetField(const vtkm::cont::DataSet& data) const;
VTKM_CONT TerminationType GetTermination(const vtkm::cont::DataSet& data) const;
VTKM_CONT AnalysisType GetAnalysis(const vtkm::cont::DataSet& data) const;
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); }
};
}
}
} // namespace vtkm::filter::flow
#endif // vtk_m_filter_flow_WarpXStreamline_h

90
vtkm/filter/flow/XGCPoincare.cxx Normal file
View File

@ -0,0 +1,90 @@
//============================================================================
// 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 <vtkm/filter/flow/XGCPoincare.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
VTKM_CONT XGCPoincare::FieldType XGCPoincare::GetField(const vtkm::cont::DataSet& dataset) const
{
const auto& field1 = this->GetActiveFieldName(0);
const auto& field2 = this->GetActiveFieldName(1);
const auto& field3 = this->GetActiveFieldName(2);
const auto& field4 = this->GetActiveFieldName(3);
const auto& field5 = this->GetActiveFieldName(4);
const auto& field6 = this->GetActiveFieldName(5);
/*
XGCField(const FieldComponentType& as_ff,
const FieldVecType& _dAs_ff_rzp,
const FieldComponentType& coeff_1D,
const FieldComponentType& coeff_2D,
const FieldComponentType& psi,
const FieldVecType& _B_rzp,
const XGCParams& params,
bool useDeltaBScale,
bool useBScale,
const Association assoc)
*/
CompArrayType as_ff;
CompArrayType coeff_1D;
CompArrayType coeff_2D;
CompArrayType psi;
VecArrayType dAs_ff_rzp;
VecArrayType b_rzp;
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(field1).GetData(), as_ff);
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(field2).GetData(), dAs_ff_rzp);
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(field3).GetData(), coeff_1D);
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(field4).GetData(), coeff_2D);
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(field5).GetData(), b_rzp);
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(field6).GetData(), psi);
vtkm::worklet::flow::XGCParams params = this->GetXGCParams();
vtkm::FloatDefault period = this->GetPeriod();
return XGCPoincare::FieldType(
as_ff, dAs_ff_rzp, coeff_1D, coeff_2D, b_rzp, psi, params, period, false, 1.0, false, 1.0);
}
VTKM_CONT XGCPoincare::TerminationType XGCPoincare::GetTermination(
const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return XGCPoincare::TerminationType(this->NumberOfSteps, this->MaxPunctures);
}
VTKM_CONT XGCPoincare::AnalysisType XGCPoincare::GetAnalysis(
const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
CompArrayType coeff_1D;
CompArrayType coeff_2D;
const auto& field3 = this->GetActiveFieldName(2);
const auto& field4 = this->GetActiveFieldName(3);
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(field3).GetData(), coeff_1D);
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(field4).GetData(), coeff_2D);
vtkm::worklet::flow::XGCParams params = this->GetXGCParams();
vtkm::FloatDefault period = this->GetPeriod();
return XGCPoincare::AnalysisType(
this->NumberOfSteps, this->MaxPunctures, coeff_1D, coeff_2D, params, period);
}
}
}
} // namespace vtkm::filter::flow

106
vtkm/filter/flow/XGCPoincare.h Normal file
View File

@ -0,0 +1,106 @@
//============================================================================
// 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_filter_flow_XGCPoincare_h
#define vtk_m_filter_flow_XGCPoincare_h
#include <vtkm/filter/flow/FilterParticleAdvectionSteadyState.h>
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
#include <vtkm/filter/flow/worklet/Analysis.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/Termination.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
class XGCPoincare;
template <>
struct FlowTraits<XGCPoincare>
{
using ParticleType = vtkm::Particle;
using TerminationType = vtkm::worklet::flow::PoincareTermination;
using AnalysisType = vtkm::worklet::flow::XGCPoincare<ParticleType>;
using CompArrayType = vtkm::cont::ArrayHandle<vtkm::FloatDefault>;
using VecArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::flow::XGCField<CompArrayType, VecArrayType>;
};
/// \brief Advect particles in a vector field.
/// Takes as input a vector field and seed locations and generates the
/// end points for each seed through the vector field.
class VTKM_FILTER_FLOW_EXPORT XGCPoincare
: public vtkm::filter::flow::FilterParticleAdvectionSteadyState<XGCPoincare>
{
public:
using ParticleType = typename FlowTraits<XGCPoincare>::ParticleType;
using TerminationType = typename FlowTraits<XGCPoincare>::TerminationType;
using AnalysisType = typename FlowTraits<XGCPoincare>::AnalysisType;
using CompArrayType = typename FlowTraits<XGCPoincare>::CompArrayType;
using VecArrayType = typename FlowTraits<XGCPoincare>::VecArrayType;
using FieldType = typename FlowTraits<XGCPoincare>::FieldType;
VTKM_CONT FieldType GetField(const vtkm::cont::DataSet& data) const;
VTKM_CONT TerminationType GetTermination(const vtkm::cont::DataSet& data) const;
VTKM_CONT AnalysisType GetAnalysis(const vtkm::cont::DataSet& data) const;
VTKM_CONT void SetXGCParams(const vtkm::worklet::flow::XGCParams& params)
{
this->Params = params;
}
VTKM_CONT vtkm::worklet::flow::XGCParams GetXGCParams() const { return this->Params; }
VTKM_CONT void SetPeriod(vtkm::FloatDefault period) { this->Period = period; }
VTKM_CONT vtkm::FloatDefault GetPeriod() const { return this->Period; }
VTKM_CONT void SetMaxPunctures(vtkm::Id maxPunctures) { this->MaxPunctures = maxPunctures; }
VTKM_CONT vtkm::Id GetMaxPunctures() const { return this->MaxPunctures; }
// As_ff
VTKM_CONT void SetField1(const std::string& name) { this->SetActiveField(0, name); }
// dAs_ff
VTKM_CONT void SetField2(const std::string& name) { this->SetActiveField(1, name); }
// Bfield
VTKM_CONT void SetField3(const std::string& name) { this->SetActiveField(2, name); }
// Psi
VTKM_CONT void SetField4(const std::string& name) { this->SetActiveField(3, name); }
// coeff_1d
VTKM_CONT void SetField5(const std::string& name) { this->SetActiveField(4, name); }
// coeff_2d
VTKM_CONT void SetField6(const std::string& name) { this->SetActiveField(5, name); }
VTKM_CONT std::string GetField1() const { return this->GetActiveFieldName(0); }
VTKM_CONT std::string GetField2() const { return this->GetActiveFieldName(1); }
VTKM_CONT std::string GetField3() const { return this->GetActiveFieldName(2); }
VTKM_CONT std::string GetField4() const { return this->GetActiveFieldName(3); }
VTKM_CONT std::string GetField5() const { return this->GetActiveFieldName(4); }
VTKM_CONT std::string GetField6() const { return this->GetActiveFieldName(5); }
private:
vtkm::worklet::flow::XGCParams Params;
vtkm::FloatDefault Period;
vtkm::Id MaxPunctures;
};
}
}
} // namespace vtkm::filter::flow
#endif // vtk_m_filter_flow_XGCPoincare_h

24
vtkm/filter/flow/internal/AdvectAlgorithm.h
View File

@ -25,10 +25,12 @@ namespace flow
namespace internal
{
template <typename DSIType, template <typename> class ResultType, typename ParticleType>
template <typename DSIType>
class AdvectAlgorithm
{
public:
using ParticleType = typename DSIType::PType;
AdvectAlgorithm(const vtkm::filter::flow::internal::BoundsMap& bm,
std::vector<DSIType>& blocks,
bool useAsyncComm)
@ -40,11 +42,8 @@ public:
{
}
void Execute(vtkm::Id maxNumSteps,
vtkm::FloatDefault stepSize,
const vtkm::cont::ArrayHandle<ParticleType>& seeds)
void Execute(const vtkm::cont::ArrayHandle<ParticleType>& seeds, vtkm::FloatDefault stepSize)
{
this->SetMaxNumberOfSteps(maxNumSteps);
this->SetStepSize(stepSize);
this->SetSeeds(seeds);
this->Go();
@ -53,19 +52,17 @@ public:
vtkm::cont::PartitionedDataSet GetOutput() const
{
vtkm::cont::PartitionedDataSet output;
for (const auto& b : this->Blocks)
{
vtkm::cont::DataSet ds;
if (b.template GetOutput<ParticleType>(ds))
if (b.GetOutput(ds))
output.AppendPartition(ds);
}
return output;
}
void SetStepSize(vtkm::FloatDefault stepSize) { this->StepSize = stepSize; }
void SetMaxNumberOfSteps(vtkm::Id numSteps) { this->MaxNumberOfSteps = numSteps; }
void SetSeeds(const vtkm::cont::ArrayHandle<ParticleType>& seeds)
{
this->ClearParticles();
@ -110,10 +107,9 @@ public:
{
//make this a pointer to avoid the copy?
auto& block = this->GetDataSet(blockId);
DSIHelperInfoType bb =
DSIHelperInfo<ParticleType>(v, this->BoundsMap, this->ParticleBlockIDsMap);
block.Advect(bb, this->StepSize, this->MaxNumberOfSteps);
numTerm = this->UpdateResult(bb.Get<DSIHelperInfo<ParticleType>>());
DSIHelperInfo<ParticleType> bb(v, this->BoundsMap, this->ParticleBlockIDsMap);
block.Advect(bb, this->StepSize);
numTerm = this->UpdateResult(bb);
}
vtkm::Id numTermMessages = 0;
@ -145,7 +141,7 @@ public:
#endif
}
DataSetIntegrator<DSIType>& GetDataSet(vtkm::Id id)
DataSetIntegrator<DSIType, ParticleType>& GetDataSet(vtkm::Id id)
{
for (auto& it : this->Blocks)
if (it.GetID() == id)

35
vtkm/filter/flow/internal/AdvectAlgorithmThreaded.h
View File

@ -28,14 +28,16 @@ namespace flow
namespace internal
{
template <typename DSIType, template <typename> class ResultType, typename ParticleType>
class AdvectAlgorithmThreaded : public AdvectAlgorithm<DSIType, ResultType, ParticleType>
template <typename DSIType>
class AdvectAlgorithmThreaded : public AdvectAlgorithm<DSIType>
{
public:
using ParticleType = typename DSIType::PType;
AdvectAlgorithmThreaded(const vtkm::filter::flow::internal::BoundsMap& bm,
std::vector<DSIType>& blocks,
bool useAsyncComm)
: AdvectAlgorithm<DSIType, ResultType, ParticleType>(bm, blocks, useAsyncComm)
: AdvectAlgorithm<DSIType>(bm, blocks)
, Done(false)
, WorkerActivate(false)
{
@ -64,8 +66,7 @@ protected:
bool GetActiveParticles(std::vector<ParticleType>& particles, vtkm::Id& blockId) override
{
std::lock_guard<std::mutex> lock(this->Mutex);
bool val = this->AdvectAlgorithm<DSIType, ResultType, ParticleType>::GetActiveParticles(
particles, blockId);
bool val = this->AdvectAlgorithm<DSIType>::GetActiveParticles(particles, blockId);
this->WorkerActivate = val;
return val;
}
@ -76,7 +77,7 @@ protected:
if (!particles.empty())
{
std::lock_guard<std::mutex> lock(this->Mutex);
this->AdvectAlgorithm<DSIType, ResultType, ParticleType>::UpdateActive(particles, idsMap);
this->AdvectAlgorithm<DSIType>::UpdateActive(particles, idsMap);
//Let workers know there is new work
this->WorkerActivateCondition.notify_all();
@ -105,7 +106,7 @@ protected:
this->WorkerActivateCondition.wait(lock, [this] { return WorkerActivate || Done; });
}
void UpdateWorkerResult(vtkm::Id blockId, DSIHelperInfoType& b)
void UpdateWorkerResult(vtkm::Id blockId, DSIHelperInfo<ParticleType>& b)
{
std::lock_guard<std::mutex> lock(this->Mutex);
auto& it = this->WorkerResults[blockId];
@ -121,9 +122,8 @@ protected:
if (this->GetActiveParticles(v, blockId))
{
auto& block = this->GetDataSet(blockId);
DSIHelperInfoType bb =
DSIHelperInfo<ParticleType>(v, this->BoundsMap, this->ParticleBlockIDsMap);
block.Advect(bb, this->StepSize, this->MaxNumberOfSteps);
DSIHelperInfo<ParticleType> bb(v, this->BoundsMap, this->ParticleBlockIDsMap);
block.Advect(bb, this->StepSize);
this->UpdateWorkerResult(blockId, bb);
}
else
@ -146,17 +146,14 @@ protected:
while (this->TotalNumTerminatedParticles < this->TotalNumParticles)
{
std::unordered_map<vtkm::Id, std::vector<DSIHelperInfoType>> workerResults;
std::unordered_map<vtkm::Id, std::vector<DSIHelperInfo<ParticleType>>> workerResults;
this->GetWorkerResults(workerResults);
vtkm::Id numTerm = 0;
for (auto& it : workerResults)
{
for (auto& r : it.second)
{
numTerm += this->UpdateResult(r.Get<DSIHelperInfo<ParticleType>>());
}
numTerm += this->UpdateResult(r);
}
vtkm::Id numTermMessages = 0;
@ -177,12 +174,12 @@ protected:
if (this->Done)
return true;
return (this->AdvectAlgorithm<DSIType, ResultType, ParticleType>::GetBlockAndWait(
syncComm, numLocalTerm) &&
return (this->AdvectAlgorithm<DSIType>::GetBlockAndWait(syncComm, numLocalTerm) &&
!this->WorkerActivate && this->WorkerResults.empty());
}
void GetWorkerResults(std::unordered_map<vtkm::Id, std::vector<DSIHelperInfoType>>& results)
void GetWorkerResults(
std::unordered_map<vtkm::Id, std::vector<DSIHelperInfo<ParticleType>>>& results)
{
results.clear();
@ -198,7 +195,7 @@ protected:
std::mutex Mutex;
bool WorkerActivate;
std::condition_variable WorkerActivateCondition;
std::unordered_map<vtkm::Id, std::vector<DSIHelperInfoType>> WorkerResults;
std::unordered_map<vtkm::Id, std::vector<DSIHelperInfo<ParticleType>>> WorkerResults;
};
}

214
vtkm/filter/flow/internal/DataSetIntegrator.h
View File

@ -102,36 +102,13 @@ public:
std::vector<vtkm::Id> TermIdx;
};
using DSIHelperInfoType =
vtkm::cont::Variant<DSIHelperInfo<vtkm::Particle>, DSIHelperInfo<vtkm::ChargedParticle>>;
template <typename Derived>
template <typename Derived, typename ParticleType>
class DataSetIntegrator
{
public:
using VelocityFieldNameType = std::string;
using ElectroMagneticFieldNameType = std::pair<std::string, std::string>;
protected:
using FieldNameType = vtkm::cont::Variant<VelocityFieldNameType, ElectroMagneticFieldNameType>;
using RType =
vtkm::cont::Variant<vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle>,
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::ChargedParticle>,
vtkm::worklet::flow::StreamlineResult<vtkm::Particle>,
vtkm::worklet::flow::StreamlineResult<vtkm::ChargedParticle>>;
public:
DataSetIntegrator(vtkm::Id id,
const FieldNameType& fieldName,
vtkm::filter::flow::IntegrationSolverType solverType,
vtkm::filter::flow::VectorFieldType vecFieldType,
vtkm::filter::flow::FlowResultType resultType)
: FieldName(fieldName)
, Id(id)
DataSetIntegrator(vtkm::Id id, vtkm::filter::flow::IntegrationSolverType solverType)
: Id(id)
, SolverType(solverType)
, VecFieldType(vecFieldType)
, AdvectionResType(resultType)
, Rank(this->Comm.rank())
{
//check that things are valid.
@ -141,57 +118,34 @@ public:
VTKM_CONT void SetCopySeedFlag(bool val) { this->CopySeedArray = val; }
VTKM_CONT
void Advect(DSIHelperInfoType& b,
vtkm::FloatDefault stepSize, //move these to member data(?)
vtkm::Id maxSteps)
void Advect(DSIHelperInfo<ParticleType>& b,
vtkm::FloatDefault stepSize) //move these to member data(?)
{
Derived* inst = static_cast<Derived*>(this);
//Cast the DSIHelperInfo<ParticleType> to the concrete type and call DoAdvect.
b.CastAndCall([&](auto& concrete) { inst->DoAdvect(concrete, stepSize, maxSteps); });
inst->DoAdvect(b, stepSize);
}
template <typename ParticleType>
VTKM_CONT bool GetOutput(vtkm::cont::DataSet& ds) const;
VTKM_CONT bool GetOutput(vtkm::cont::DataSet& dataset) const
{
Derived* inst = static_cast<Derived*>(this);
return inst->GetOutput(dataset);
}
protected:
template <typename ParticleType, template <typename> class ResultType>
VTKM_CONT void UpdateResult(const ResultType<ParticleType>& result,
DSIHelperInfo<ParticleType>& dsiInfo);
VTKM_CONT bool IsParticleAdvectionResult() const
{
return this->AdvectionResType == FlowResultType::PARTICLE_ADVECT_TYPE;
}
VTKM_CONT bool IsStreamlineResult() const
{
return this->AdvectionResType == FlowResultType::STREAMLINE_TYPE;
}
template <typename ParticleType>
VTKM_CONT inline void ClassifyParticles(const vtkm::cont::ArrayHandle<ParticleType>& particles,
DSIHelperInfo<ParticleType>& dsiInfo) const;
//Data members.
vtkm::cont::Variant<VelocityFieldNameType, ElectroMagneticFieldNameType> FieldName;
vtkm::Id Id;
vtkm::filter::flow::IntegrationSolverType SolverType;
vtkm::filter::flow::VectorFieldType VecFieldType;
vtkm::filter::flow::FlowResultType AdvectionResType =
vtkm::filter::flow::FlowResultType::UNKNOWN_TYPE;
vtkmdiy::mpi::communicator Comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
vtkm::Id Rank;
bool CopySeedArray = false;
std::vector<RType> Results;
};
template <typename Derived>
template <typename ParticleType>
VTKM_CONT inline void DataSetIntegrator<Derived>::ClassifyParticles(
//template <typename Derived, typename ParticleType, typename AnalysisType>
template <typename Derived, typename ParticleType>
VTKM_CONT inline void DataSetIntegrator<Derived, ParticleType>::ClassifyParticles(
const vtkm::cont::ArrayHandle<ParticleType>& particles,
DSIHelperInfo<ParticleType>& dsiInfo) const
{
@ -275,146 +229,6 @@ VTKM_CONT inline void DataSetIntegrator<Derived>::ClassifyParticles(
dsiInfo.Validate(n);
}
template <typename Derived>
template <typename ParticleType, template <typename> class ResultType>
VTKM_CONT inline void DataSetIntegrator<Derived>::UpdateResult(
const ResultType<ParticleType>& result,
DSIHelperInfo<ParticleType>& dsiInfo)
{
this->ClassifyParticles(result.Particles, dsiInfo);
if (this->IsParticleAdvectionResult())
{
if (dsiInfo.TermIdx.empty())
return;
using ResType = vtkm::worklet::flow::ParticleAdvectionResult<ParticleType>;
auto indicesAH = vtkm::cont::make_ArrayHandle(dsiInfo.TermIdx, vtkm::CopyFlag::Off);
auto termPerm = vtkm::cont::make_ArrayHandlePermutation(indicesAH, result.Particles);
vtkm::cont::ArrayHandle<ParticleType> termParticles;
vtkm::cont::Algorithm::Copy(termPerm, termParticles);
ResType termRes(termParticles);
this->Results.emplace_back(termRes);
}
else if (this->IsStreamlineResult())
this->Results.emplace_back(result);
}
template <typename Derived>
template <typename ParticleType>
VTKM_CONT inline bool DataSetIntegrator<Derived>::GetOutput(vtkm::cont::DataSet& ds) const
{
std::size_t nResults = this->Results.size();
if (nResults == 0)
return false;
if (this->IsParticleAdvectionResult())
{
using ResType = vtkm::worklet::flow::ParticleAdvectionResult<ParticleType>;
std::vector<vtkm::cont::ArrayHandle<ParticleType>> allParticles;
allParticles.reserve(nResults);
for (const auto& vres : this->Results)
allParticles.emplace_back(vres.template Get<ResType>().Particles);
vtkm::cont::ArrayHandle<vtkm::Vec3f> pts;
vtkm::cont::ParticleArrayCopy(allParticles, pts);
vtkm::Id numPoints = pts.GetNumberOfValues();
if (numPoints > 0)
{
//Create coordinate system and vertex cell set.
ds.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", pts));
vtkm::cont::CellSetSingleType<> cells;
vtkm::cont::ArrayHandleIndex conn(numPoints);
vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
vtkm::cont::ArrayCopy(conn, connectivity);
cells.Fill(numPoints, vtkm::CELL_SHAPE_VERTEX, 1, connectivity);
ds.SetCellSet(cells);
}
}
else if (this->IsStreamlineResult())
{
using ResType = vtkm::worklet::flow::StreamlineResult<ParticleType>;
//Easy case with one result.
if (nResults == 1)
{
const auto& res = this->Results[0].template Get<ResType>();
ds.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", res.Positions));
ds.SetCellSet(res.PolyLines);
}
else
{
std::vector<vtkm::Id> posOffsets(nResults, 0);
vtkm::Id totalNumCells = 0, totalNumPts = 0;
for (std::size_t i = 0; i < nResults; i++)
{
const auto& res = this->Results[i].template Get<ResType>();
if (i == 0)
posOffsets[i] = 0;
else
posOffsets[i] = totalNumPts;
totalNumPts += res.Positions.GetNumberOfValues();
totalNumCells += res.PolyLines.GetNumberOfCells();
}
//Append all the points together.
vtkm::cont::ArrayHandle<vtkm::Vec3f> appendPts;
appendPts.Allocate(totalNumPts);
for (std::size_t i = 0; i < nResults; i++)
{
const auto& res = this->Results[i].template Get<ResType>();
// copy all values into appendPts starting at offset.
vtkm::cont::Algorithm::CopySubRange(
res.Positions, 0, res.Positions.GetNumberOfValues(), appendPts, posOffsets[i]);
}
ds.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", appendPts));
//Create polylines.
std::vector<vtkm::Id> numPtsPerCell(static_cast<std::size_t>(totalNumCells));
std::size_t off = 0;
for (std::size_t i = 0; i < nResults; i++)
{
const auto& res = this->Results[i].template Get<ResType>();
vtkm::Id nCells = res.PolyLines.GetNumberOfCells();
for (vtkm::Id j = 0; j < nCells; j++)
numPtsPerCell[off++] = static_cast<vtkm::Id>(res.PolyLines.GetNumberOfPointsInCell(j));
}
auto numPointsPerCellArray = vtkm::cont::make_ArrayHandle(numPtsPerCell, vtkm::CopyFlag::Off);
vtkm::cont::ArrayHandle<vtkm::Id> cellIndex;
vtkm::Id connectivityLen =
vtkm::cont::Algorithm::ScanExclusive(numPointsPerCellArray, cellIndex);
vtkm::cont::ArrayHandleIndex connCount(connectivityLen);
vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
vtkm::cont::ArrayCopy(connCount, connectivity);
vtkm::cont::ArrayHandle<vtkm::UInt8> cellTypes;
auto polyLineShape = vtkm::cont::make_ArrayHandleConstant<vtkm::UInt8>(
vtkm::CELL_SHAPE_POLY_LINE, totalNumCells);
vtkm::cont::ArrayCopy(polyLineShape, cellTypes);
auto offsets = vtkm::cont::ConvertNumComponentsToOffsets(numPointsPerCellArray);
vtkm::cont::CellSetExplicit<> polyLines;
polyLines.Fill(totalNumPts, cellTypes, connectivity, offsets);
ds.SetCellSet(polyLines);
}
}
else
{
throw vtkm::cont::ErrorFilterExecution("Unsupported ParticleAdvectionResultType");
}
return true;
}
}
}
}

340
vtkm/filter/flow/internal/DataSetIntegratorSteadyState.h
View File

@ -22,252 +22,150 @@ namespace flow
{
namespace internal
{
class DataSetIntegratorSteadyState
: public vtkm::filter::flow::internal::DataSetIntegrator<DataSetIntegratorSteadyState>
{
public:
DataSetIntegratorSteadyState(const vtkm::cont::DataSet& ds,
vtkm::Id id,
const FieldNameType& fieldName,
vtkm::filter::flow::IntegrationSolverType solverType,
vtkm::filter::flow::VectorFieldType vecFieldType,
vtkm::filter::flow::FlowResultType resultType)
: vtkm::filter::flow::internal::DataSetIntegrator<DataSetIntegratorSteadyState>(id,
fieldName,
solverType,
vecFieldType,
resultType)
, DataSet(ds)
{
}
VTKM_CONT inline void DoAdvect(DSIHelperInfo<vtkm::Particle>& b,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps);
VTKM_CONT inline void DoAdvect(DSIHelperInfo<vtkm::ChargedParticle>& b,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps);
protected:
template <typename ArrayType>
VTKM_CONT void GetVelocityField(
vtkm::worklet::flow::VelocityField<ArrayType>& velocityField) const
{
if (this->FieldName.IsType<VelocityFieldNameType>())
{
const auto& fieldNm = this->FieldName.Get<VelocityFieldNameType>();
auto assoc = this->DataSet.GetField(fieldNm).GetAssociation();
ArrayType arr;
vtkm::cont::ArrayCopyShallowIfPossible(this->DataSet.GetField(fieldNm).GetData(), arr);
velocityField = vtkm::worklet::flow::VelocityField<ArrayType>(arr, assoc);
}
else
throw vtkm::cont::ErrorFilterExecution("Velocity field vector type not available");
}
template <typename ArrayType>
VTKM_CONT void GetElectroMagneticField(
vtkm::worklet::flow::ElectroMagneticField<ArrayType>& ebField) const
{
if (this->FieldName.IsType<ElectroMagneticFieldNameType>())
{
const auto& fieldNm = this->FieldName.Get<ElectroMagneticFieldNameType>();
const auto& electric = fieldNm.first;
const auto& magnetic = fieldNm.second;
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<ArrayType>(eField, bField, eAssoc);
}
else
throw vtkm::cont::ErrorFilterExecution("Electromagnetic field vector type not available");
}
private:
vtkm::cont::DataSet DataSet;
};
namespace internal
{
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using VelocityFieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
using VelocityEvalType = vtkm::worklet::flow::GridEvaluator<VelocityFieldType>;
using EBFieldType = vtkm::worklet::flow::ElectroMagneticField<ArrayType>;
using EBEvalType = vtkm::worklet::flow::GridEvaluator<EBFieldType>;
//template <typename GridEvalType, typename ParticleType>
//class AdvectHelper;
template <typename FieldType, typename ParticleType>
class AdvectHelper //<FieldType, ParticleType>
namespace detail
{
template <typename ParticleType,
typename FieldType,
typename TerminationType,
typename AnalysisType>
class AdvectHelperSteadyState
{
public:
using WorkletType = vtkm::worklet::flow::ParticleAdvection;
using SteadyStateGridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
static void Advect(const FieldType& vecField,
const vtkm::cont::DataSet& ds,
vtkm::cont::ArrayHandle<ParticleType>& seedArray,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
const IntegrationSolverType& solverType,
vtkm::worklet::flow::ParticleAdvectionResult<ParticleType>& result)
{
if (solverType == IntegrationSolverType::RK4_TYPE)
{
DoAdvect<vtkm::worklet::flow::ParticleAdvection,
vtkm::worklet::flow::ParticleAdvectionResult,
vtkm::worklet::flow::RK4Integrator>(
vecField, ds, seedArray, stepSize, maxSteps, result);
}
else if (solverType == IntegrationSolverType::EULER_TYPE)
{
DoAdvect<vtkm::worklet::flow::ParticleAdvection,
vtkm::worklet::flow::ParticleAdvectionResult,
vtkm::worklet::flow::EulerIntegrator>(
vecField, ds, seedArray, stepSize, maxSteps, result);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported Integrator type");
}
static void Advect(const FieldType& vecField,
const vtkm::cont::DataSet& ds,
vtkm::cont::ArrayHandle<ParticleType>& seedArray,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
const IntegrationSolverType& solverType,
vtkm::worklet::flow::StreamlineResult<ParticleType>& result)
{
if (solverType == IntegrationSolverType::RK4_TYPE)
{
DoAdvect<vtkm::worklet::flow::Streamline,
vtkm::worklet::flow::StreamlineResult,
vtkm::worklet::flow::RK4Integrator>(
vecField, ds, seedArray, stepSize, maxSteps, result);
}
else if (solverType == IntegrationSolverType::EULER_TYPE)
{
DoAdvect<vtkm::worklet::flow::Streamline,
vtkm::worklet::flow::StreamlineResult,
vtkm::worklet::flow::EulerIntegrator>(
vecField, ds, seedArray, stepSize, maxSteps, result);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported Integrator type");
}
template <typename WorkletType,
template <typename>
class ResultType,
template <typename>
class SolverType>
static void DoAdvect(const FieldType& vecField,
const vtkm::cont::DataSet& ds,
vtkm::cont::ArrayHandle<ParticleType>& seedArray,
template <template <typename> class SolverType>
static void DoAdvect(vtkm::cont::ArrayHandle<ParticleType>& seedArray,
const FieldType& field,
const vtkm::cont::DataSet& dataset,
const TerminationType& termination,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
ResultType<ParticleType>& result)
AnalysisType& analysis)
{
using StepperType =
vtkm::worklet::flow::Stepper<SolverType<SteadyStateGridEvalType>, SteadyStateGridEvalType>;
SteadyStateGridEvalType eval(dataset, field);
StepperType stepper(eval, stepSize);
WorkletType worklet;
SteadyStateGridEvalType eval(ds, vecField);
StepperType stepper(eval, stepSize);
result = worklet.Run(stepper, seedArray, maxSteps);
worklet.Run(stepper, seedArray, termination, analysis);
}
static void Advect(vtkm::cont::ArrayHandle<ParticleType>& seedArray,
const FieldType& field,
const vtkm::cont::DataSet& dataset,
const TerminationType& termination,
const IntegrationSolverType& solverType,
vtkm::FloatDefault stepSize,
AnalysisType& analysis)
{
if (solverType == IntegrationSolverType::RK4_TYPE)
{
DoAdvect<vtkm::worklet::flow::RK4Integrator>(
seedArray, field, dataset, termination, stepSize, analysis);
}
else if (solverType == IntegrationSolverType::EULER_TYPE)
{
DoAdvect<vtkm::worklet::flow::EulerIntegrator>(
seedArray, field, dataset, termination, stepSize, analysis);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported Integrator type");
}
};
}
} // namespace detail
VTKM_CONT inline void DataSetIntegratorSteadyState::DoAdvect(DSIHelperInfo<vtkm::Particle>& b,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps)
template <typename ParticleType,
typename FieldType,
typename TerminationType,
typename AnalysisType>
class DataSetIntegratorSteadyState
: public vtkm::filter::flow::internal::DataSetIntegrator<
DataSetIntegratorSteadyState<ParticleType, FieldType, TerminationType, AnalysisType>,
ParticleType>
{
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
public:
using BaseType = vtkm::filter::flow::internal::DataSetIntegrator<
DataSetIntegratorSteadyState<ParticleType, FieldType, TerminationType, AnalysisType>,
ParticleType>;
using PType = ParticleType;
using FType = FieldType;
using TType = TerminationType;
using AType = AnalysisType;
auto copyFlag = (this->CopySeedArray ? vtkm::CopyFlag::On : vtkm::CopyFlag::Off);
auto seedArray = vtkm::cont::make_ArrayHandle(b.Particles, copyFlag);
if (this->VecFieldType == VectorFieldType::VELOCITY_FIELD_TYPE)
DataSetIntegratorSteadyState(vtkm::Id id,
const FieldType& field,
const vtkm::cont::DataSet& dataset,
vtkm::filter::flow::IntegrationSolverType solverType,
const TerminationType& termination,
const AnalysisType& analysis)
: BaseType(id, solverType)
, Field(field)
, Dataset(dataset)
, Termination(termination)
, Analysis(analysis)
{
using FieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
FieldType vecField;
this->GetVelocityField(vecField);
}
using AHType = internal::AdvectHelper<internal::VelocityFieldType, vtkm::Particle>;
VTKM_CONT inline void DoAdvect(vtkm::filter::flow::internal::DSIHelperInfo<ParticleType>& block,
vtkm::FloatDefault stepSize)
{
auto copyFlag = (this->CopySeedArray ? vtkm::CopyFlag::On : vtkm::CopyFlag::Off);
auto seedArray = vtkm::cont::make_ArrayHandle(block.V, copyFlag);
if (this->IsParticleAdvectionResult())
using AdvectionHelper =
detail::AdvectHelperSteadyState<ParticleType, FieldType, TerminationType, AnalysisType>;
AnalysisType analysis;
analysis.UseAsTemplate(this->Analysis);
AdvectionHelper::Advect(seedArray,
this->Field,
this->Dataset,
this->Termination,
this->SolverType,
stepSize,
analysis);
this->UpdateResult(analysis, block);
}
VTKM_CONT void UpdateResult(AnalysisType& analysis,
vtkm::filter::flow::internal::DSIHelperInfo<ParticleType>& dsiInfo)
{
this->ClassifyParticles(analysis.Particles, dsiInfo);
if (std::is_same<AnalysisType, vtkm::worklet::flow::NoAnalysis<ParticleType>>::value)
{
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> result;
AHType::Advect(
vecField, this->DataSet, seedArray, stepSize, maxSteps, this->SolverType, result);
this->UpdateResult(result, b);
}
else if (this->IsStreamlineResult())
{
vtkm::worklet::flow::StreamlineResult<vtkm::Particle> result;
AHType::Advect(
vecField, this->DataSet, seedArray, stepSize, maxSteps, this->SolverType, result);
this->UpdateResult(result, b);
if (dsiInfo.TermIdx.empty())
return;
auto indicesAH = vtkm::cont::make_ArrayHandle(dsiInfo.TermIdx, vtkm::CopyFlag::Off);
auto termPerm = vtkm::cont::make_ArrayHandlePermutation(indicesAH, analysis.Particles);
vtkm::cont::ArrayHandle<ParticleType> termParticles;
vtkm::cont::Algorithm::Copy(termPerm, termParticles);
analysis.FinalizeAnalysis(termParticles);
this->Analyses.emplace_back(analysis);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported result type");
{
this->Analyses.emplace_back(analysis);
}
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported vector field type");
}
VTKM_CONT inline void DataSetIntegratorSteadyState::DoAdvect(
DSIHelperInfo<vtkm::ChargedParticle>& b,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps)
{
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
auto copyFlag = (this->CopySeedArray ? vtkm::CopyFlag::On : vtkm::CopyFlag::Off);
auto seedArray = vtkm::cont::make_ArrayHandle(b.Particles, copyFlag);
if (this->VecFieldType == VectorFieldType::ELECTRO_MAGNETIC_FIELD_TYPE)
VTKM_CONT bool GetOutput(vtkm::cont::DataSet& ds) const
{
using FieldType = vtkm::worklet::flow::ElectroMagneticField<ArrayType>;
FieldType ebField;
this->GetElectroMagneticField(ebField);
using AHType = internal::AdvectHelper<internal::EBFieldType, vtkm::ChargedParticle>;
if (this->IsParticleAdvectionResult())
{
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::ChargedParticle> result;
AHType::Advect(
ebField, this->DataSet, seedArray, stepSize, maxSteps, this->SolverType, result);
this->UpdateResult(result, b);
}
else if (this->IsStreamlineResult())
{
vtkm::worklet::flow::StreamlineResult<vtkm::ChargedParticle> result;
AHType::Advect(
ebField, this->DataSet, seedArray, stepSize, maxSteps, this->SolverType, result);
this->UpdateResult(result, b);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported result type");
std::size_t nAnalyses = this->Analyses.size();
if (nAnalyses == 0)
return false;
return AnalysisType::MakeDataSet(ds, this->Analyses);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported vector field type");
}
private:
FieldType Field;
vtkm::cont::DataSet Dataset;
TerminationType Termination;
// Used as a template to initialize successive analysis objects.
AnalysisType Analysis;
std::vector<AnalysisType> Analyses;
};
}
}

351
vtkm/filter/flow/internal/DataSetIntegratorUnsteadyState.h
View File

@ -23,242 +23,173 @@ namespace flow
{
namespace internal
{
class DataSetIntegratorUnsteadyState
: public vtkm::filter::flow::internal::DataSetIntegrator<DataSetIntegratorUnsteadyState>
namespace detail
{
template <typename ParticleType,
typename FieldType,
typename TerminationType,
typename AnalysisType>
class AdvectHelperUnsteadyState
{
public:
DataSetIntegratorUnsteadyState(const vtkm::cont::DataSet& ds1,
using WorkletType = vtkm::worklet::flow::ParticleAdvection;
using UnsteadyStateGridEvalType = vtkm::worklet::flow::TemporalGridEvaluator<FieldType>;
template <template <typename> class SolverType>
static void DoAdvect(vtkm::cont::ArrayHandle<ParticleType>& seedArray,
const FieldType& field1,
const vtkm::cont::DataSet& ds1,
vtkm::FloatDefault t1,
const FieldType& field2,
const vtkm::cont::DataSet& ds2,
vtkm::FloatDefault t2,
const TerminationType& termination,
vtkm::FloatDefault stepSize,
AnalysisType& analysis)
{
using StepperType = vtkm::worklet::flow::Stepper<SolverType<UnsteadyStateGridEvalType>,
UnsteadyStateGridEvalType>;
WorkletType worklet;
UnsteadyStateGridEvalType eval(ds1, t1, field1, ds2, t2, field2);
StepperType stepper(eval, stepSize);
worklet.Run(stepper, seedArray, termination, analysis);
}
static void Advect(vtkm::cont::ArrayHandle<ParticleType>& seedArray,
const FieldType& field1,
const vtkm::cont::DataSet& ds1,
vtkm::FloatDefault t1,
const FieldType& field2,
const vtkm::cont::DataSet& ds2,
vtkm::FloatDefault t2,
const TerminationType& termination,
const IntegrationSolverType& solverType,
vtkm::FloatDefault stepSize,
AnalysisType& analysis)
{
if (solverType == IntegrationSolverType::RK4_TYPE)
{
DoAdvect<vtkm::worklet::flow::RK4Integrator>(
seedArray, field1, ds1, t1, field2, ds2, t2, termination, stepSize, analysis);
}
else if (solverType == IntegrationSolverType::EULER_TYPE)
{
DoAdvect<vtkm::worklet::flow::EulerIntegrator>(
seedArray, field1, ds1, t1, field2, ds2, t2, termination, stepSize, analysis);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported Integrator type");
}
};
} //namespace detail
template <typename ParticleType,
typename FieldType,
typename TerminationType,
typename AnalysisType>
class DataSetIntegratorUnsteadyState
: public vtkm::filter::flow::internal::DataSetIntegrator<
DataSetIntegratorUnsteadyState<ParticleType, FieldType, TerminationType, AnalysisType>,
ParticleType>
{
public:
using BaseType = vtkm::filter::flow::internal::DataSetIntegrator<
DataSetIntegratorUnsteadyState<ParticleType, FieldType, TerminationType, AnalysisType>,
ParticleType>;
using PType = ParticleType;
using FType = FieldType;
using TType = TerminationType;
using AType = AnalysisType;
DataSetIntegratorUnsteadyState(vtkm::Id id,
const FieldType& field1,
const FieldType& field2,
const vtkm::cont::DataSet& ds1,
const vtkm::cont::DataSet& ds2,
vtkm::FloatDefault t1,
vtkm::FloatDefault t2,
vtkm::Id id,
const vtkm::filter::flow::internal::DataSetIntegrator<
DataSetIntegratorUnsteadyState>::FieldNameType& fieldName,
vtkm::filter::flow::IntegrationSolverType solverType,
vtkm::filter::flow::VectorFieldType vecFieldType,
vtkm::filter::flow::FlowResultType resultType)
: vtkm::filter::flow::internal::DataSetIntegrator<DataSetIntegratorUnsteadyState>(id,
fieldName,
solverType,
vecFieldType,
resultType)
const TerminationType& termination,
const AnalysisType& analysis)
: BaseType(id, solverType)
, Field1(field1)
, Field2(field2)
, DataSet1(ds1)
, DataSet2(ds2)
, Time1(t1)
, Time2(t2)
, Termination(termination)
, Analysis(analysis)
{
}
VTKM_CONT inline void DoAdvect(DSIHelperInfo<vtkm::Particle>& b,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps);
VTKM_CONT inline void DoAdvect(DSIHelperInfo<vtkm::ChargedParticle>& b,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps);
protected:
template <typename ArrayType>
VTKM_CONT void GetVelocityFields(
vtkm::worklet::flow::VelocityField<ArrayType>& velocityField1,
vtkm::worklet::flow::VelocityField<ArrayType>& velocityField2) const
VTKM_CONT inline void DoAdvect(vtkm::filter::flow::internal::DSIHelperInfo<ParticleType>& block,
vtkm::FloatDefault stepSize)
{
if (this->FieldName.GetIndex() == this->FieldName.GetIndexOf<VelocityFieldNameType>())
auto copyFlag = (this->CopySeedArray ? vtkm::CopyFlag::On : vtkm::CopyFlag::Off);
auto seedArray = vtkm::cont::make_ArrayHandle(block.V, copyFlag);
using AdvectionHelper =
detail::AdvectHelperUnsteadyState<ParticleType, FieldType, TerminationType, AnalysisType>;
AnalysisType analysis;
analysis.UseAsTemplate(this->Analysis);
AdvectionHelper::Advect(seedArray,
this->Field1,
this->DataSet1,
this->Time1,
this->Field2,
this->DataSet2,
this->Time2,
this->Termination,
this->SolverType,
stepSize,
analysis);
this->UpdateResult(analysis, block);
}
VTKM_CONT void UpdateResult(AnalysisType& analysis,
vtkm::filter::flow::internal::DSIHelperInfo<ParticleType>& dsiInfo)
{
this->ClassifyParticles(analysis.Particles, dsiInfo);
if (std::is_same<AnalysisType, vtkm::worklet::flow::NoAnalysis<ParticleType>>::value)
{
const auto& fieldNm = this->FieldName.Get<VelocityFieldNameType>();
auto assoc = this->DataSet1.GetField(fieldNm).GetAssociation();
if (assoc != this->DataSet2.GetField(fieldNm).GetAssociation())
throw vtkm::cont::ErrorFilterExecution(
"Unsteady state velocity fields have differnt associations");
ArrayType arr1, arr2;
vtkm::cont::ArrayCopyShallowIfPossible(this->DataSet1.GetField(fieldNm).GetData(), arr1);
vtkm::cont::ArrayCopyShallowIfPossible(this->DataSet2.GetField(fieldNm).GetData(), arr2);
velocityField1 = vtkm::worklet::flow::VelocityField<ArrayType>(arr1, assoc);
velocityField2 = vtkm::worklet::flow::VelocityField<ArrayType>(arr2, assoc);
if (dsiInfo.TermIdx.empty())
return;
auto indicesAH = vtkm::cont::make_ArrayHandle(dsiInfo.TermIdx, vtkm::CopyFlag::Off);
auto termPerm = vtkm::cont::make_ArrayHandlePermutation(indicesAH, analysis.Particles);
vtkm::cont::ArrayHandle<ParticleType> termParticles;
vtkm::cont::Algorithm::Copy(termPerm, termParticles);
analysis.FinalizeAnalysis(termParticles);
this->Analyses.emplace_back(analysis);
}
else
throw vtkm::cont::ErrorFilterExecution("Velocity field vector type not available");
{
this->Analyses.emplace_back(analysis);
}
}
VTKM_CONT bool GetOutput(vtkm::cont::DataSet& ds) const
{
std::size_t nAnalyses = this->Analyses.size();
if (nAnalyses == 0)
return false;
return AnalysisType::MakeDataSet(ds, this->Analyses);
}
private:
FieldType Field1;
FieldType Field2;
vtkm::cont::DataSet DataSet1;
vtkm::cont::DataSet DataSet2;
vtkm::FloatDefault Time1;
vtkm::FloatDefault Time2;
TerminationType Termination;
AnalysisType Analysis;
std::vector<AnalysisType> Analyses;
};
namespace internal
{
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using VelocityFieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
using UnsteadyStateGridEvalType = vtkm::worklet::flow::TemporalGridEvaluator<VelocityFieldType>;
template <typename GridEvalType, typename ParticleType>
class AdvectHelper;
template <typename ParticleType>
class AdvectHelper<UnsteadyStateGridEvalType, ParticleType>
{
public:
static void Advect(const VelocityFieldType& velField1,
const vtkm::cont::DataSet& ds1,
vtkm::FloatDefault t1,
const VelocityFieldType& velField2,
const vtkm::cont::DataSet& ds2,
vtkm::FloatDefault t2,
vtkm::cont::ArrayHandle<ParticleType>& seedArray,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
const IntegrationSolverType& solverType,
vtkm::worklet::flow::ParticleAdvectionResult<ParticleType>& result)
{
if (solverType == IntegrationSolverType::RK4_TYPE)
{
DoAdvect<vtkm::worklet::flow::ParticleAdvection,
vtkm::worklet::flow::ParticleAdvectionResult,
vtkm::worklet::flow::RK4Integrator>(
velField1, ds1, t1, velField2, ds2, t2, seedArray, stepSize, maxSteps, result);
}
else if (solverType == IntegrationSolverType::EULER_TYPE)
{
DoAdvect<vtkm::worklet::flow::ParticleAdvection,
vtkm::worklet::flow::ParticleAdvectionResult,
vtkm::worklet::flow::EulerIntegrator>(
velField1, ds1, t1, velField2, ds2, t2, seedArray, stepSize, maxSteps, result);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported Integrator type");
}
static void Advect(const VelocityFieldType& velField1,
const vtkm::cont::DataSet& ds1,
vtkm::FloatDefault t1,
const VelocityFieldType& velField2,
const vtkm::cont::DataSet& ds2,
vtkm::FloatDefault t2,
vtkm::cont::ArrayHandle<ParticleType>& seedArray,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
const IntegrationSolverType& solverType,
vtkm::worklet::flow::StreamlineResult<ParticleType>& result)
{
if (solverType == IntegrationSolverType::RK4_TYPE)
{
DoAdvect<vtkm::worklet::flow::Streamline,
vtkm::worklet::flow::StreamlineResult,
vtkm::worklet::flow::RK4Integrator>(
velField1, ds1, t1, velField2, ds2, t2, seedArray, stepSize, maxSteps, result);
}
else if (solverType == IntegrationSolverType::EULER_TYPE)
{
DoAdvect<vtkm::worklet::flow::Streamline,
vtkm::worklet::flow::StreamlineResult,
vtkm::worklet::flow::EulerIntegrator>(
velField1, ds1, t1, velField2, ds2, t2, seedArray, stepSize, maxSteps, result);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported Integrator type");
}
template <typename WorkletType,
template <typename>
class ResultType,
template <typename>
class SolverType>
static void DoAdvect(const VelocityFieldType& velField1,
const vtkm::cont::DataSet& ds1,
vtkm::FloatDefault t1,
const VelocityFieldType& velField2,
const vtkm::cont::DataSet& ds2,
vtkm::FloatDefault t2,
vtkm::cont::ArrayHandle<ParticleType>& seedArray,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
ResultType<ParticleType>& result)
{
using StepperType = vtkm::worklet::flow::Stepper<SolverType<UnsteadyStateGridEvalType>,
UnsteadyStateGridEvalType>;
WorkletType worklet;
UnsteadyStateGridEvalType eval(ds1, t1, velField1, ds2, t2, velField2);
StepperType stepper(eval, stepSize);
result = worklet.Run(stepper, seedArray, maxSteps);
}
};
}
VTKM_CONT inline void DataSetIntegratorUnsteadyState::DoAdvect(DSIHelperInfo<vtkm::Particle>& b,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps)
{
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
auto copyFlag = (this->CopySeedArray ? vtkm::CopyFlag::On : vtkm::CopyFlag::Off);
auto seedArray = vtkm::cont::make_ArrayHandle(b.Particles, copyFlag);
using AHType = internal::AdvectHelper<internal::UnsteadyStateGridEvalType, vtkm::Particle>;
if (this->VecFieldType == VectorFieldType::VELOCITY_FIELD_TYPE)
{
using FieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
FieldType velField1, velField2;
this->GetVelocityFields(velField1, velField2);
if (this->IsParticleAdvectionResult())
{
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> result;
AHType::Advect(velField1,
this->DataSet1,
this->Time1,
velField2,
this->DataSet2,
this->Time2,
seedArray,
stepSize,
maxSteps,
this->SolverType,
result);
this->UpdateResult(result, b);
}
else if (this->IsStreamlineResult())
{
vtkm::worklet::flow::StreamlineResult<vtkm::Particle> result;
AHType::Advect(velField1,
this->DataSet1,
this->Time1,
velField2,
this->DataSet2,
this->Time2,
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");
}
VTKM_CONT inline void DataSetIntegratorUnsteadyState::DoAdvect(
DSIHelperInfo<vtkm::ChargedParticle>& vtkmNotUsed(b),
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");
}
}
}
}

84
vtkm/filter/flow/internal/ParticleAdvector.h
View File

@ -29,93 +29,47 @@ template <typename DSIType>
class ParticleAdvector
{
public:
using ParticleType = typename DSIType::PType;
ParticleAdvector(const vtkm::filter::flow::internal::BoundsMap& bm,
const std::vector<DSIType>& blocks,
const bool& useThreaded,
const bool& useAsyncComm,
const vtkm::filter::flow::FlowResultType& parType)
const bool& useAsyncComm)
: Blocks(blocks)
, BoundsMap(bm)
, ResultType(parType)
, UseAsynchronousCommunication(useAsyncComm)
, UseThreadedAlgorithm(useThreaded)
{
}
vtkm::cont::PartitionedDataSet Execute(vtkm::Id numSteps,
vtkm::FloatDefault stepSize,
const vtkm::cont::UnknownArrayHandle& seeds)
{
using ParticleTypeList = vtkm::List<vtkm::Particle, vtkm::ChargedParticle>;
vtkm::cont::PartitionedDataSet result;
seeds.CastAndCallForTypes<ParticleTypeList, VTKM_DEFAULT_STORAGE_LIST>(
[&](const auto& concreteSeeds) {
result = this->Execute(numSteps, stepSize, concreteSeeds);
});
return result;
}
private:
template <typename AlgorithmType, typename ParticleType>
vtkm::cont::PartitionedDataSet RunAlgo(vtkm::Id numSteps,
vtkm::FloatDefault stepSize,
const vtkm::cont::ArrayHandle<ParticleType>& seeds)
{
AlgorithmType algo(this->BoundsMap, this->Blocks, this->UseAsynchronousCommunication);
algo.Execute(numSteps, stepSize, seeds);
return algo.GetOutput();
}
template <typename ParticleType>
vtkm::cont::PartitionedDataSet Execute(vtkm::Id numSteps,
vtkm::FloatDefault stepSize,
const vtkm::cont::ArrayHandle<ParticleType>& seeds)
vtkm::cont::PartitionedDataSet Execute(const vtkm::cont::ArrayHandle<ParticleType>& seeds,
vtkm::FloatDefault stepSize)
{
if (!this->UseThreadedAlgorithm)
{
if (this->ResultType == vtkm::filter::flow::FlowResultType::PARTICLE_ADVECT_TYPE)
{
using AlgorithmType = vtkm::filter::flow::internal::
AdvectAlgorithm<DSIType, vtkm::worklet::flow::ParticleAdvectionResult, ParticleType>;
return this->RunAlgo<AlgorithmType, ParticleType>(numSteps, stepSize, seeds);
}
else
{
using AlgorithmType = vtkm::filter::flow::internal::
AdvectAlgorithm<DSIType, vtkm::worklet::flow::StreamlineResult, ParticleType>;
return this->RunAlgo<AlgorithmType, ParticleType>(numSteps, stepSize, seeds);
}
using AlgorithmType = vtkm::filter::flow::internal::AdvectAlgorithm<DSIType>;
return this->RunAlgo<AlgorithmType>(seeds, stepSize);
}
else
{
if (this->ResultType == vtkm::filter::flow::FlowResultType::PARTICLE_ADVECT_TYPE)
{
using AlgorithmType = vtkm::filter::flow::internal::AdvectAlgorithmThreaded<
DSIType,
vtkm::worklet::flow::ParticleAdvectionResult,
ParticleType>;
return this->RunAlgo<AlgorithmType, ParticleType>(numSteps, stepSize, seeds);
}
else
{
using AlgorithmType = vtkm::filter::flow::internal::
AdvectAlgorithmThreaded<DSIType, vtkm::worklet::flow::StreamlineResult, ParticleType>;
return this->RunAlgo<AlgorithmType, ParticleType>(numSteps, stepSize, seeds);
}
using AlgorithmType = vtkm::filter::flow::internal::AdvectAlgorithmThreaded<DSIType>;
return this->RunAlgo<AlgorithmType>(seeds, stepSize);
}
}
private:
template <typename AlgorithmType>
vtkm::cont::PartitionedDataSet RunAlgo(const vtkm::cont::ArrayHandle<ParticleType>& seeds,
vtkm::FloatDefault stepSize)
{
AlgorithmType algo(this->BoundsMap, this->Blocks);
algo.Execute(seeds, stepSize);
return algo.GetOutput();
}
std::vector<DSIType> Blocks;
vtkm::filter::flow::internal::BoundsMap BoundsMap;
FlowResultType ResultType;
bool UseAsynchronousCommunication = true;
std::vector<DSIType> Blocks;
bool UseThreadedAlgorithm;
};

2
vtkm/filter/flow/testing/CMakeLists.txt
View File

@ -28,7 +28,7 @@ if (TARGET vtkm_rendering_testing)
endif()
vtkm_unit_tests(
SOURCES ${filter_unit_tests}
DEVICE_SOURCES ${filter_unit_tests}
DEVICE_SOURCES ${worklet_unit_tests}
USE_VTKM_JOB_POOL
)

4
vtkm/filter/flow/testing/UnitTestStreamlineFilter.cxx
View File

@ -372,7 +372,6 @@ void TestPartitionedDataSet(vtkm::Id num, bool useGhost, FilterType fType)
seedArray = vtkm::cont::make_ArrayHandle({ vtkm::Particle(vtkm::Vec3f(.2f, 1.0f, .2f), 0),
vtkm::Particle(vtkm::Vec3f(.2f, 2.0f, .2f), 1) });
vtkm::Id numSeeds = seedArray.GetNumberOfValues();
if (fType == FilterType::STREAMLINE || fType == FilterType::PATHLINE)
{
vtkm::cont::PartitionedDataSet out;
@ -382,7 +381,6 @@ void TestPartitionedDataSet(vtkm::Id num, bool useGhost, FilterType fType)
streamline.SetStepSize(0.1f);
streamline.SetNumberOfSteps(100000);
streamline.SetSeeds(seedArray);
streamline.SetActiveField(fieldName);
out = streamline.Execute(pds);
}
@ -590,8 +588,6 @@ void TestStreamlineFilters()
FilterType::PATHLINE,
FilterType::PATH_PARTICLE };
//fTypes = {FilterType::PARTICLE_ADVECTION,FilterType::STREAMLINE};
//fTypes = {FilterType::PATHLINE};
for (int n = 1; n < 3; n++)
{
for (auto useGhost : flags)

10
vtkm/filter/flow/testing/UnitTestStreamlineFilterWarpX.cxx
View File

@ -11,7 +11,7 @@
#include <vtkm/Particle.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/flow/Streamline.h>
#include <vtkm/filter/flow/WarpXStreamline.h>
#include <vtkm/io/VTKDataSetReader.h>
namespace
@ -43,7 +43,7 @@ void GenerateChargedParticles(const vtkm::cont::ArrayHandle<vtkm::Vec3f>& pos,
}
}
void TestStreamlineFilters()
void TestFilters()
{
std::string particleFile = vtkm::cont::testing::Testing::DataPath("misc/warpXparticles.vtk");
std::string fieldFile = vtkm::cont::testing::Testing::DataPath("misc/warpXfields.vtk");
@ -88,12 +88,10 @@ void TestStreamlineFilters()
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;
vtkm::filter::flow::WarpXStreamline 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");
@ -109,5 +107,5 @@ void TestStreamlineFilters()
int UnitTestStreamlineFilterWarpX(int argc, char* argv[])
{
return vtkm::cont::testing::Testing::Run(TestStreamlineFilters, argc, argv);
return vtkm::cont::testing::Testing::Run(TestFilters, argc, argv);
}

125
vtkm/filter/flow/testing/UnitTestWorkletParticleAdvection.cxx
View File

@ -14,6 +14,7 @@
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/flow/worklet/Analysis.h>
#include <vtkm/filter/flow/worklet/EulerIntegrator.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/GridEvaluators.h>
@ -21,6 +22,7 @@
#include <vtkm/filter/flow/worklet/Particles.h>
#include <vtkm/filter/flow/worklet/RK4Integrator.h>
#include <vtkm/filter/flow/worklet/Stepper.h>
#include <vtkm/filter/flow/worklet/Termination.h>
#include <vtkm/filter/mesh_info/GhostCellClassify.h>
#include <vtkm/io/VTKDataSetReader.h>
#include <vtkm/worklet/testing/GenerateTestDataSets.h>
@ -369,6 +371,8 @@ void TestGhostCellEvaluators()
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
using Termination = vtkm::worklet::flow::NormalTermination;
using Analysis = vtkm::worklet::flow::NoAnalysis<vtkm::Particle>;
constexpr vtkm::Id nX = 6;
constexpr vtkm::Id nY = 6;
@ -407,9 +411,11 @@ void TestGhostCellEvaluators()
seeds.push_back(vtkm::Particle(vtkm::Vec3f(3, 3, 3), 3));
auto seedArray = vtkm::cont::make_ArrayHandle(seeds, vtkm::CopyFlag::Off);
auto res = pa.Run(rk4, seedArray, 10000);
Termination termination(10000);
Analysis analysis;
pa.Run(rk4, seedArray, termination, analysis);
auto posPortal = res.Particles.ReadPortal();
auto posPortal = analysis.Particles.ReadPortal();
vtkm::Id numSeeds = seedArray.GetNumberOfValues();
for (vtkm::Id i = 0; i < numSeeds; i++)
{
@ -428,10 +434,9 @@ void TestGhostCellEvaluators()
}
}
void ValidateParticleAdvectionResult(
const vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle>& res,
vtkm::Id nSeeds,
vtkm::Id maxSteps)
void ValidateParticleAdvectionResult(const vtkm::worklet::flow::NoAnalysis<vtkm::Particle>& res,
vtkm::Id nSeeds,
vtkm::Id maxSteps)
{
VTKM_TEST_ASSERT(res.Particles.GetNumberOfValues() == nSeeds,
"Number of output particles does not match input.");
@ -449,7 +454,7 @@ void ValidateParticleAdvectionResult(
}
}
void ValidateStreamlineResult(const vtkm::worklet::flow::StreamlineResult<vtkm::Particle>& res,
void ValidateStreamlineResult(const vtkm::worklet::flow::StreamlineAnalysis<vtkm::Particle>& res,
vtkm::Id nSeeds,
vtkm::Id maxSteps)
{
@ -471,6 +476,9 @@ void TestIntegrators()
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::flow::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using Termination = vtkm::worklet::flow::NormalTermination;
using Analysis = vtkm::worklet::flow::NoAnalysis<vtkm::Particle>;
const vtkm::Id3 dims(5, 5, 5);
const vtkm::Bounds bounds(0., 1., 0., 1., .0, .1);
@ -496,22 +504,25 @@ void TestIntegrators()
GenerateRandomParticles(points, 3, bounds);
vtkm::worklet::flow::ParticleAdvection pa;
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> res;
Termination termination(maxSteps);
Analysis analysis;
{
auto seeds = vtkm::cont::make_ArrayHandle(points, vtkm::CopyFlag::On);
using IntegratorType = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<IntegratorType, GridEvalType>;
Stepper rk4(eval, stepSize);
res = pa.Run(rk4, seeds, maxSteps);
ValidateParticleAdvectionResult(res, nSeeds, maxSteps);
pa.Run(rk4, seeds, termination, analysis);
ValidateParticleAdvectionResult(analysis, nSeeds, maxSteps);
}
{
auto seeds = vtkm::cont::make_ArrayHandle(points, vtkm::CopyFlag::On);
using IntegratorType = vtkm::worklet::flow::EulerIntegrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<IntegratorType, GridEvalType>;
Stepper euler(eval, stepSize);
res = pa.Run(euler, seeds, maxSteps);
ValidateParticleAdvectionResult(res, nSeeds, maxSteps);
pa.Run(euler, seeds, termination, analysis);
ValidateParticleAdvectionResult(analysis, nSeeds, maxSteps);
//res = pa.Run(euler, seeds, maxSteps);
//ValidateParticleAdvectionResult(res, nSeeds, maxSteps);
}
}
}
@ -523,6 +534,9 @@ void TestParticleWorkletsWithDataSetTypes()
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
using Termination = vtkm::worklet::flow::NormalTermination;
using PAnalysis = vtkm::worklet::flow::NoAnalysis<vtkm::Particle>;
using SAnalysis = vtkm::worklet::flow::StreamlineAnalysis<vtkm::Particle>;
vtkm::FloatDefault stepSize = 0.01f;
const vtkm::Id3 dims(5, 5, 5);
@ -574,34 +588,36 @@ void TestParticleWorkletsWithDataSetTypes()
if (i < 2)
{
vtkm::worklet::flow::ParticleAdvection pa;
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> res;
Termination termination(maxSteps);
PAnalysis analysis;
if (i == 0)
{
auto seeds = vtkm::cont::make_ArrayHandle(pts, vtkm::CopyFlag::On);
res = pa.Run(rk4, seeds, maxSteps);
pa.Run(rk4, seeds, termination, analysis);
}
else
{
auto seeds = vtkm::cont::make_ArrayHandle(pts2, vtkm::CopyFlag::On);
res = pa.Run(rk4, seeds, maxSteps);
pa.Run(rk4, seeds, termination, analysis);
}
ValidateParticleAdvectionResult(res, nSeeds, maxSteps);
ValidateParticleAdvectionResult(analysis, nSeeds, maxSteps);
}
else
{
vtkm::worklet::flow::Streamline s;
vtkm::worklet::flow::StreamlineResult<vtkm::Particle> res;
vtkm::worklet::flow::ParticleAdvection pa;
Termination termination(maxSteps);
SAnalysis analysis(maxSteps);
if (i == 2)
{
auto seeds = vtkm::cont::make_ArrayHandle(pts, vtkm::CopyFlag::On);
res = s.Run(rk4, seeds, maxSteps);
pa.Run(rk4, seeds, termination, analysis);
}
else
{
auto seeds = vtkm::cont::make_ArrayHandle(pts2, vtkm::CopyFlag::On);
res = s.Run(rk4, seeds, maxSteps);
pa.Run(rk4, seeds, termination, analysis);
}
ValidateStreamlineResult(res, nSeeds, maxSteps);
ValidateStreamlineResult(analysis, nSeeds, maxSteps);
}
}
}
@ -629,6 +645,8 @@ void TestParticleStatus()
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
using Termination = vtkm::worklet::flow::NormalTermination;
using Analysis = vtkm::worklet::flow::NoAnalysis<vtkm::Particle>;
vtkm::Id maxSteps = 1000;
vtkm::FloatDefault stepSize = 0.01f;
@ -644,7 +662,10 @@ void TestParticleStatus()
pts.push_back(vtkm::Particle(vtkm::Vec3f(-1, -1, -1), 1));
auto seedsArray = vtkm::cont::make_ArrayHandle(pts, vtkm::CopyFlag::On);
pa.Run(rk4, seedsArray, maxSteps);
Termination termination(maxSteps);
Analysis analysis;
pa.Run(rk4, seedsArray, termination, analysis);
auto portal = seedsArray.ReadPortal();
bool tookStep0 = portal.Get(0).GetStatus().CheckTookAnySteps();
@ -679,6 +700,9 @@ void TestWorkletsBasic()
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
using Termination = vtkm::worklet::flow::NormalTermination;
using PAnalysis = vtkm::worklet::flow::NoAnalysis<vtkm::Particle>;
using SAnalysis = vtkm::worklet::flow::StreamlineAnalysis<vtkm::Particle>;
vtkm::FloatDefault stepSize = 0.01f;
const vtkm::Id3 dims(5, 5, 5);
@ -738,15 +762,15 @@ void TestWorkletsBasic()
if (w == "particleAdvection")
{
vtkm::worklet::flow::ParticleAdvection pa;
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> res;
res = pa.Run(rk4, seedsArray, maxSteps);
Termination termination(maxSteps);
PAnalysis analysis;
pa.Run(rk4, seedsArray, termination, analysis);
vtkm::Id numRequiredPoints = static_cast<vtkm::Id>(endPts.size());
VTKM_TEST_ASSERT(res.Particles.GetNumberOfValues() == numRequiredPoints,
VTKM_TEST_ASSERT(analysis.Particles.GetNumberOfValues() == numRequiredPoints,
"Wrong number of points in particle advection result.");
auto portal = res.Particles.ReadPortal();
for (vtkm::Id i = 0; i < res.Particles.GetNumberOfValues(); i++)
auto portal = analysis.Particles.ReadPortal();
for (vtkm::Id i = 0; i < analysis.Particles.GetNumberOfValues(); i++)
{
VTKM_TEST_ASSERT(portal.Get(i).GetPosition() == endPts[static_cast<std::size_t>(i)],
"Particle advection point is wrong");
@ -762,18 +786,19 @@ void TestWorkletsBasic()
}
else if (w == "streamline")
{
vtkm::worklet::flow::Streamline s;
vtkm::worklet::flow::StreamlineResult<vtkm::Particle> res;
res = s.Run(rk4, seedsArray, maxSteps);
vtkm::worklet::flow::ParticleAdvection pa;
Termination termination(maxSteps);
SAnalysis analysis(maxSteps);
pa.Run(rk4, seedsArray, termination, analysis);
vtkm::Id numRequiredPoints = static_cast<vtkm::Id>(samplePts.size());
VTKM_TEST_ASSERT(res.Positions.GetNumberOfValues() == numRequiredPoints,
VTKM_TEST_ASSERT(analysis.Streams.GetNumberOfValues() == numRequiredPoints,
"Wrong number of points in streamline result.");
//Make sure all the points match.
auto parPortal = res.Particles.ReadPortal();
for (vtkm::Id i = 0; i < res.Particles.GetNumberOfValues(); i++)
auto parPortal = analysis.Particles.ReadPortal();
for (vtkm::Id i = 0; i < analysis.Particles.GetNumberOfValues(); i++)
{
VTKM_TEST_ASSERT(parPortal.Get(i).GetPosition() == endPts[static_cast<std::size_t>(i)],
"Streamline end point is wrong");
@ -786,19 +811,19 @@ void TestWorkletsBasic()
"Streamline particle did not terminate");
}
auto posPortal = res.Positions.ReadPortal();
for (vtkm::Id i = 0; i < res.Positions.GetNumberOfValues(); i++)
auto posPortal = analysis.Streams.ReadPortal();
for (vtkm::Id i = 0; i < analysis.Streams.GetNumberOfValues(); i++)
VTKM_TEST_ASSERT(posPortal.Get(i) == samplePts[static_cast<std::size_t>(i)],
"Streamline points do not match");
vtkm::Id numCells = res.PolyLines.GetNumberOfCells();
vtkm::Id numCells = analysis.PolyLines.GetNumberOfCells();
VTKM_TEST_ASSERT(numCells == static_cast<vtkm::Id>(pts.size()),
"Wrong number of polylines in streamline");
for (vtkm::Id i = 0; i < numCells; i++)
{
VTKM_TEST_ASSERT(res.PolyLines.GetCellShape(i) == vtkm::CELL_SHAPE_POLY_LINE,
VTKM_TEST_ASSERT(analysis.PolyLines.GetCellShape(i) == vtkm::CELL_SHAPE_POLY_LINE,
"Wrong cell type in streamline.");
VTKM_TEST_ASSERT(res.PolyLines.GetNumberOfPointsInCell(i) ==
VTKM_TEST_ASSERT(analysis.PolyLines.GetNumberOfPointsInCell(i) ==
static_cast<vtkm::Id>(maxSteps + 1),
"Wrong number of points in streamline cell");
}
@ -869,6 +894,9 @@ void TestParticleAdvectionFile(const std::string& fileName,
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
using Termination = vtkm::worklet::flow::NormalTermination;
using PAnalysis = vtkm::worklet::flow::NoAnalysis<vtkm::Particle>;
using SAnalysis = vtkm::worklet::flow::StreamlineAnalysis<vtkm::Particle>;
VTKM_TEST_ASSERT(ds.HasField(fieldName), "Data set missing a field named ", fieldName);
vtkm::cont::Field& field = ds.GetField(fieldName);
@ -887,6 +915,7 @@ void TestParticleAdvectionFile(const std::string& fileName,
FieldType velocities(fieldArray);
GridEvalType eval(ds.GetCoordinateSystem(), ds.GetCellSet(), velocities);
Stepper rk4(eval, stepSize);
Termination termination(maxSteps);
for (int i = 0; i < 2; i++)
{
@ -898,18 +927,16 @@ void TestParticleAdvectionFile(const std::string& fileName,
if (i == 0)
{
vtkm::worklet::flow::ParticleAdvection pa;
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> res;
res = pa.Run(rk4, seedArray, maxSteps);
ValidateResult(res, maxSteps, endPts);
PAnalysis analysis;
pa.Run(rk4, seedArray, termination, analysis);
ValidateResult(analysis, maxSteps, endPts);
}
else if (i == 1)
{
vtkm::worklet::flow::Streamline s;
vtkm::worklet::flow::StreamlineResult<vtkm::Particle> res;
res = s.Run(rk4, seedArray, maxSteps);
ValidateResult(res, maxSteps, endPts);
vtkm::worklet::flow::ParticleAdvection pa;
SAnalysis analysis(maxSteps);
pa.Run(rk4, seedArray, termination, analysis);
ValidateResult(analysis, maxSteps, endPts);
}
}
}

11
vtkm/filter/flow/worklet/Analysis.cxx Normal file
View File

@ -0,0 +1,11 @@
//============================================================================
// 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 <vtkm/filter/flow/worklet/Analysis.h>

412
vtkm/filter/flow/worklet/Analysis.h Normal file
View File

@ -0,0 +1,412 @@
//=============================================================================
//
// 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 vtkm_worklet_particleadvection_analysis
#define vtkm_worklet_particleadvection_analysis
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ConvertNumComponentsToOffsets.h>
#include <vtkm/cont/Invoker.h>
#include <vtkm/cont/ParticleArrayCopy.h>
#include <vtkm/worklet/WorkletMapField.h>
namespace vtkm
{
namespace worklet
{
namespace flow
{
template <typename ParticleType>
class NoAnalysisExec
{
public:
VTKM_EXEC_CONT
NoAnalysisExec() {}
//template <typename ParticleType>
VTKM_EXEC void Analyze(const vtkm::Id index,
const ParticleType& oldParticle,
const ParticleType& newParticle)
{
// Do nothing
(void)index;
(void)oldParticle;
(void)newParticle;
}
};
template <typename ParticleType>
class NoAnalysis : public vtkm::cont::ExecutionObjectBase
{
public:
// Intended to store advected particles after Finalize
vtkm::cont::ArrayHandle<ParticleType> Particles;
VTKM_CONT
NoAnalysis()
: Particles()
{
}
VTKM_CONT
void UseAsTemplate(const NoAnalysis& other) { (void)other; }
VTKM_CONT
//template <typename ParticleType>
void InitializeAnalysis(const vtkm::cont::ArrayHandle<ParticleType>& particles)
{
(void)particles;
}
VTKM_CONT
//template <typename ParticleType>
void FinalizeAnalysis(vtkm::cont::ArrayHandle<ParticleType>& particles)
{
this->Particles = particles; //, vtkm::CopyFlag::Off);
}
VTKM_CONT NoAnalysisExec<ParticleType> PrepareForExecution(vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
{
(void)device;
(void)token;
return NoAnalysisExec<ParticleType>();
}
VTKM_CONT bool SupportPushOutOfBounds() const { return true; }
VTKM_CONT static bool MakeDataSet(vtkm::cont::DataSet& dataset,
const std::vector<NoAnalysis>& results)
{
size_t nResults = results.size();
std::vector<vtkm::cont::ArrayHandle<ParticleType>> allParticles;
allParticles.reserve(nResults);
for (const auto& vres : results)
allParticles.emplace_back(vres.Particles);
vtkm::cont::ArrayHandle<vtkm::Vec3f> pts;
vtkm::cont::ParticleArrayCopy(allParticles, pts);
vtkm::Id numPoints = pts.GetNumberOfValues();
if (numPoints > 0)
{
//Create coordinate system and vertex cell set.
dataset.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", pts));
vtkm::cont::CellSetSingleType<> cells;
vtkm::cont::ArrayHandleIndex conn(numPoints);
vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
vtkm::cont::ArrayCopy(conn, connectivity);
cells.Fill(numPoints, vtkm::CELL_SHAPE_VERTEX, 1, connectivity);
dataset.SetCellSet(cells);
}
return true;
}
};
template <typename ParticleType>
class StreamlineAnalysisExec
{
public:
VTKM_EXEC_CONT
StreamlineAnalysisExec()
: NumParticles(0)
, MaxSteps(0)
, Streams()
, StreamLengths()
, Validity()
{
}
VTKM_CONT
StreamlineAnalysisExec(vtkm::Id numParticles,
vtkm::Id maxSteps,
const vtkm::cont::ArrayHandle<vtkm::Vec3f>& streams,
const vtkm::cont::ArrayHandle<vtkm::Id>& streamLengths,
const vtkm::cont::ArrayHandle<vtkm::Id>& validity,
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token)
: NumParticles(numParticles)
, MaxSteps(maxSteps + 1)
{
Streams = streams.PrepareForOutput(this->NumParticles * this->MaxSteps, device, token);
StreamLengths = streamLengths.PrepareForInPlace(device, token);
Validity = validity.PrepareForInPlace(device, token);
}
//template <typename ParticleType>
VTKM_EXEC void Analyze(const vtkm::Id index,
const ParticleType& oldParticle,
const ParticleType& newParticle)
{
vtkm::Id loc = index * MaxSteps;
vtkm::Id streamLength = this->StreamLengths.Get(index);
if (streamLength == 0)
{
this->StreamLengths.Set(index, 1);
this->Streams.Set(loc, oldParticle.GetPosition());
this->Validity.Set(loc, 1);
++streamLength;
}
loc += streamLength;
this->StreamLengths.Set(index, ++streamLength);
this->Streams.Set(loc, newParticle.GetPosition());
this->Validity.Set(loc, 1);
}
private:
using IdPortal = typename vtkm::cont::ArrayHandle<vtkm::Id>::WritePortalType;
using VecPortal = typename vtkm::cont::ArrayHandle<vtkm::Vec3f>::WritePortalType;
vtkm::Id NumParticles;
vtkm::Id MaxSteps;
VecPortal Streams;
IdPortal StreamLengths;
IdPortal Validity;
};
namespace detail
{
class GetSteps : public vtkm::worklet::WorkletMapField
{
public:
VTKM_CONT
GetSteps() {}
using ControlSignature = void(FieldIn, FieldOut);
using ExecutionSignature = void(_1, _2);
template <typename ParticleType>
VTKM_EXEC void operator()(const ParticleType& p, vtkm::Id& numSteps) const
{
numSteps = p.GetNumberOfSteps();
}
};
class ComputeNumPoints : public vtkm::worklet::WorkletMapField
{
public:
VTKM_CONT
ComputeNumPoints() {}
using ControlSignature = void(FieldIn, FieldIn, FieldOut);
using ExecutionSignature = void(_1, _2, _3);
// Offset is number of points in streamline.
// 1 (inital point) + number of steps taken (p.NumSteps - initalNumSteps)
template <typename ParticleType>
VTKM_EXEC void operator()(const ParticleType& p,
const vtkm::Id& initialNumSteps,
vtkm::Id& diff) const
{
diff = 1 + p.GetNumberOfSteps() - initialNumSteps;
}
};
} // namespace detail
template <typename ParticleType>
class StreamlineAnalysis : public vtkm::cont::ExecutionObjectBase
{
public:
// Intended to store advected particles after Finalize
vtkm::cont::ArrayHandle<ParticleType> Particles;
vtkm::cont::ArrayHandle<vtkm::Vec3f> Streams;
vtkm::cont::CellSetExplicit<> PolyLines;
//Helper functor for compacting history
struct IsOne
{
template <typename T>
VTKM_EXEC_CONT bool operator()(const T& x) const
{
return x == T(1);
}
};
VTKM_CONT
StreamlineAnalysis()
: Particles()
, MaxSteps(0)
{
}
VTKM_CONT
StreamlineAnalysis(vtkm::Id maxSteps)
: Particles()
, MaxSteps(maxSteps)
{
}
VTKM_CONT
void UseAsTemplate(const StreamlineAnalysis& other) { this->MaxSteps = other.MaxSteps; }
VTKM_CONT StreamlineAnalysisExec<ParticleType> PrepareForExecution(
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
{
return StreamlineAnalysisExec<ParticleType>(this->NumParticles,
this->MaxSteps,
this->Streams,
this->StreamLengths,
this->Validity,
device,
token);
}
VTKM_CONT bool SupportPushOutOfBounds() const { return true; }
VTKM_CONT
//template <typename ParticleType>
void InitializeAnalysis(const vtkm::cont::ArrayHandle<ParticleType>& particles)
{
this->NumParticles = particles.GetNumberOfValues();
//Create ValidPointArray initialized to zero.
vtkm::cont::ArrayHandleConstant<vtkm::Id> validity(0,
this->NumParticles * (this->MaxSteps + 1));
vtkm::cont::ArrayCopy(validity, this->Validity);
//Create StepCountArray initialized to zero.
vtkm::cont::ArrayHandleConstant<vtkm::Id> streamLengths(0, this->NumParticles);
vtkm::cont::ArrayCopy(streamLengths, this->StreamLengths);
// Initialize InitLengths
vtkm::Id numSeeds = static_cast<vtkm::Id>(particles.GetNumberOfValues());
vtkm::cont::ArrayHandleIndex idxArray(numSeeds);
vtkm::cont::Invoker invoker;
invoker(detail::GetSteps{}, particles, this->InitialLengths);
}
VTKM_CONT
//template <typename ParticleType>
void FinalizeAnalysis(vtkm::cont::ArrayHandle<ParticleType>& particles)
{
vtkm::Id numSeeds = particles.GetNumberOfValues();
vtkm::cont::ArrayHandle<vtkm::Vec3f> positions;
vtkm::cont::Algorithm::CopyIf(this->Streams, this->Validity, positions, IsOne());
vtkm::cont::Algorithm::Copy(positions, this->Streams);
// Create the cells
vtkm::cont::ArrayHandle<vtkm::Id> numPoints;
vtkm::cont::Invoker invoker;
invoker(detail::ComputeNumPoints{}, particles, this->InitialLengths, numPoints);
vtkm::cont::ArrayHandle<vtkm::Id> cellIndex;
vtkm::Id connectivityLen = vtkm::cont::Algorithm::ScanExclusive(numPoints, cellIndex);
vtkm::cont::ArrayHandleIndex connCount(connectivityLen);
vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
vtkm::cont::ArrayCopy(connCount, connectivity);
vtkm::cont::ArrayHandle<vtkm::UInt8> cellTypes;
auto polyLineShape =
vtkm::cont::make_ArrayHandleConstant<vtkm::UInt8>(vtkm::CELL_SHAPE_POLY_LINE, numSeeds);
vtkm::cont::ArrayCopy(polyLineShape, cellTypes);
auto offsets = vtkm::cont::ConvertNumComponentsToOffsets(numPoints);
this->PolyLines.Fill(this->Streams.GetNumberOfValues(), cellTypes, connectivity, offsets);
this->Particles = particles;
//this->Validity.ReleaseResources();
//this->InitialLengths.ReleaseResources();
//this->StreamLengths.ReleaseResources();
}
VTKM_CONT static bool MakeDataSet(vtkm::cont::DataSet& dataset,
const std::vector<StreamlineAnalysis>& results)
{
size_t nResults = results.size();
if (nResults == 1)
{
const auto& res = results[0];
dataset.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", res.Streams));
dataset.SetCellSet(res.PolyLines);
}
else
{
std::vector<vtkm::Id> posOffsets(nResults, 0);
vtkm::Id totalNumCells = 0, totalNumPts = 0;
for (std::size_t i = 0; i < nResults; i++)
{
const auto& res = results[i];
if (i == 0)
posOffsets[i] = 0;
else
posOffsets[i] = totalNumPts;
totalNumPts += res.Streams.GetNumberOfValues();
totalNumCells += res.PolyLines.GetNumberOfCells();
}
//Append all the points together.
vtkm::cont::ArrayHandle<vtkm::Vec3f> appendPts;
appendPts.Allocate(totalNumPts);
for (std::size_t i = 0; i < nResults; i++)
{
const auto& res = results[i];
// copy all values into appendPts starting at offset.
vtkm::cont::Algorithm::CopySubRange(
res.Streams, 0, res.Streams.GetNumberOfValues(), appendPts, posOffsets[i]);
}
dataset.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", appendPts));
//Create polylines.
std::vector<vtkm::Id> numPtsPerCell(static_cast<std::size_t>(totalNumCells));
std::size_t off = 0;
for (std::size_t i = 0; i < nResults; i++)
{
const auto& res = results[i];
vtkm::Id nCells = res.PolyLines.GetNumberOfCells();
for (vtkm::Id j = 0; j < nCells; j++)
numPtsPerCell[off++] = static_cast<vtkm::Id>(res.PolyLines.GetNumberOfPointsInCell(j));
}
auto numPointsPerCellArray = vtkm::cont::make_ArrayHandle(numPtsPerCell, vtkm::CopyFlag::Off);
vtkm::cont::ArrayHandle<vtkm::Id> cellIndex;
vtkm::Id connectivityLen =
vtkm::cont::Algorithm::ScanExclusive(numPointsPerCellArray, cellIndex);
vtkm::cont::ArrayHandleIndex connCount(connectivityLen);
vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
vtkm::cont::ArrayCopy(connCount, connectivity);
vtkm::cont::ArrayHandle<vtkm::UInt8> cellTypes;
auto polyLineShape = vtkm::cont::make_ArrayHandleConstant<vtkm::UInt8>(
vtkm::CELL_SHAPE_POLY_LINE, totalNumCells);
vtkm::cont::ArrayCopy(polyLineShape, cellTypes);
auto offsets = vtkm::cont::ConvertNumComponentsToOffsets(numPointsPerCellArray);
vtkm::cont::CellSetExplicit<> polyLines;
polyLines.Fill(totalNumPts, cellTypes, connectivity, offsets);
dataset.SetCellSet(polyLines);
}
return true;
}
private:
vtkm::Id NumParticles;
vtkm::Id MaxSteps;
vtkm::cont::ArrayHandle<vtkm::Id> StreamLengths;
vtkm::cont::ArrayHandle<vtkm::Id> InitialLengths;
vtkm::cont::ArrayHandle<vtkm::Id> Validity;
};
} // namespace particleadvection
} // namespace worklet
} // namespace vtkm
#include <vtkm/filter/flow/worklet/XGCPoincare.h>
#endif //vtkm_worklet_particleadvection_analysis

8
vtkm/filter/flow/worklet/CMakeLists.txt
View File

@ -9,7 +9,7 @@
##============================================================================
set(headers
ParticleAdvection.h
Analysis.h
CellInterpolationHelper.h
EulerIntegrator.h
Field.h
@ -18,11 +18,17 @@ set(headers
IntegratorStatus.h
LagrangianStructures.h
Particles.h
ParticleAdvection.h
ParticleAdvectionWorklets.h
RK4Integrator.h
TemporalGridEvaluators.h
Stepper.h
StreamSurface.h
TemporalGridEvaluators.h
Termination.h
XGCHelper.h
XGCField.h
XGCPoincare.h
)
#-----------------------------------------------------------------------------

2
vtkm/filter/flow/worklet/EulerIntegrator.h
View File

@ -34,7 +34,7 @@ public:
}
template <typename Particle>
VTKM_EXEC IntegratorStatus CheckStep(Particle& particle,
VTKM_EXEC IntegratorStatus CheckStep(const Particle& particle,
vtkm::FloatDefault stepLength,
vtkm::Vec3f& velocity) const
{

102
vtkm/filter/flow/worklet/Field.h
View File

@ -31,6 +31,7 @@ class ExecutionVelocityField
public:
using FieldPortalType = typename FieldArrayType::ReadPortalType;
using Association = vtkm::cont::Field::Association;
using DelegateToField = std::false_type;
VTKM_CONT
ExecutionVelocityField(const FieldArrayType& velocityValues,
@ -68,17 +69,67 @@ public:
value = vtkm::make_Vec(velocityInterp);
}
template <typename Point, typename Locator, typename Helper>
VTKM_EXEC bool GetValue(const Point& vtkmNotUsed(point),
const vtkm::FloatDefault& vtkmNotUsed(time),
vtkm::VecVariable<Point, 2>& vtkmNotUsed(out),
const Locator& vtkmNotUsed(locator),
const Helper& vtkmNotUsed(helper)) const
{
//TODO Raise Error : Velocity Field should not allow this path
return false;
}
private:
FieldPortalType VelocityValues;
Association Assoc;
};
template <typename FieldArrayType>
class VelocityField : public vtkm::cont::ExecutionObjectBase
{
public:
using ExecutionType = ExecutionVelocityField<FieldArrayType>;
using Association = vtkm::cont::Field::Association;
VTKM_CONT
VelocityField() = default;
VTKM_CONT
VelocityField(const FieldArrayType& fieldValues)
: FieldValues(fieldValues)
, Assoc(vtkm::cont::Field::Association::Points)
{
}
VTKM_CONT
VelocityField(const FieldArrayType& fieldValues, const Association assoc)
: FieldValues(fieldValues)
, Assoc(assoc)
{
if (assoc != Association::Points && assoc != Association::Cells)
throw("Unsupported field association");
}
VTKM_CONT
const ExecutionType PrepareForExecution(vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
{
return ExecutionType(this->FieldValues, this->Assoc, device, token);
}
private:
FieldArrayType FieldValues;
Association Assoc;
};
template <typename FieldArrayType>
class ExecutionElectroMagneticField
{
public:
using FieldPortalType = typename FieldArrayType::ReadPortalType;
using Association = vtkm::cont::Field::Association;
using DelegateToField = std::false_type;
VTKM_CONT
ExecutionElectroMagneticField(const FieldArrayType& electricValues,
@ -124,50 +175,23 @@ public:
value = vtkm::make_Vec(electricInterp, magneticInterp);
}
template <typename Point, typename Locator, typename Helper>
VTKM_EXEC bool GetValue(const Point& vtkmNotUsed(point),
const vtkm::FloatDefault& vtkmNotUsed(time),
vtkm::VecVariable<Point, 2>& vtkmNotUsed(out),
const Locator& vtkmNotUsed(locator),
const Helper& vtkmNotUsed(helper)) const
{
//TODO : Raise Error : Velocity Field should not allow this path
return false;
}
private:
FieldPortalType ElectricValues;
FieldPortalType MagneticValues;
Association Assoc;
};
template <typename FieldArrayType>
class VelocityField : public vtkm::cont::ExecutionObjectBase
{
public:
using ExecutionType = ExecutionVelocityField<FieldArrayType>;
using Association = vtkm::cont::Field::Association;
VTKM_CONT
VelocityField() = default;
VTKM_CONT
VelocityField(const FieldArrayType& fieldValues)
: FieldValues(fieldValues)
, Assoc(vtkm::cont::Field::Association::Points)
{
}
VTKM_CONT
VelocityField(const FieldArrayType& fieldValues, const Association assoc)
: FieldValues(fieldValues)
, Assoc(assoc)
{
if (assoc != Association::Points && assoc != Association::Cells)
throw("Unsupported field association");
}
VTKM_CONT
const ExecutionType PrepareForExecution(vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
{
return ExecutionType(this->FieldValues, this->Assoc, device, token);
}
private:
FieldArrayType FieldValues;
Association Assoc;
};
template <typename FieldArrayType>
class ElectroMagneticField : public vtkm::cont::ExecutionObjectBase
{
@ -213,4 +237,6 @@ private:
}
} //vtkm::worklet::flow
#include <vtkm/filter/flow/worklet/XGCField.h>
#endif //vtkm_filter_flow_worklet_Field_h

41
vtkm/filter/flow/worklet/GridEvaluators.h
View File

@ -35,6 +35,7 @@ template <typename FieldType>
class ExecutionGridEvaluator
{
using GhostCellArrayType = vtkm::cont::ArrayHandle<vtkm::UInt8>;
using ExecFieldType = typename FieldType::ExecutionType;
public:
VTKM_CONT
@ -86,9 +87,9 @@ public:
}
template <typename Point>
VTKM_EXEC GridEvaluatorStatus Evaluate(const Point& point,
const vtkm::FloatDefault& time,
vtkm::VecVariable<Point, 2>& out) const
VTKM_EXEC GridEvaluatorStatus HelpEvaluate(const Point& point,
const vtkm::FloatDefault& time,
vtkm::VecVariable<Point, 2>& out) const
{
vtkm::Id cellId = -1;
Point parametric;
@ -134,10 +135,40 @@ public:
status.SetOk();
}
return status;
}
template <typename Point>
VTKM_EXEC GridEvaluatorStatus DeligateEvaluateToField(const Point& point,
const vtkm::FloatDefault& time,
vtkm::VecVariable<Point, 2>& out) const
{
GridEvaluatorStatus status;
status.SetOk();
// TODO: Allow for getting status from deligated work from Field
if (!this->Field.GetValue(point, time, out, this->Locator, this->InterpolationHelper))
{
status.SetFail();
status.SetSpatialBounds();
}
return status;
}
template <typename Point>
VTKM_EXEC GridEvaluatorStatus Evaluate(const Point& point,
const vtkm::FloatDefault& time,
vtkm::VecVariable<Point, 2>& out) const
{
if (!ExecFieldType::DelegateToField::value)
{
return this->HelpEvaluate(point, time, out);
}
else
{
return this->DeligateEvaluateToField(point, time, out);
}
}
private:
VTKM_EXEC bool InGhostCell(const vtkm::Id& cellId) const
{
@ -150,7 +181,7 @@ private:
using GhostCellPortal = typename vtkm::cont::ArrayHandle<vtkm::UInt8>::ReadPortalType;
vtkm::Bounds Bounds;
typename FieldType::ExecutionType Field;
ExecFieldType Field;
GhostCellPortal GhostCells;
bool HaveGhostCells;
vtkm::exec::CellInterpolationHelper InterpolationHelper;

57
vtkm/filter/flow/worklet/ParticleAdvection.h
View File

@ -47,7 +47,7 @@ public:
};
} //detail
/*
template <typename ParticleType>
struct ParticleAdvectionResult
{
@ -63,12 +63,12 @@ struct ParticleAdvectionResult
vtkm::cont::ArrayHandle<ParticleType> Particles;
};
*/
class ParticleAdvection
{
public:
ParticleAdvection() {}
/*
template <typename IntegratorType, typename ParticleType, typename ParticleStorage>
void Run2(const IntegratorType& it,
vtkm::cont::ArrayHandle<ParticleType, ParticleStorage>& particles,
@ -80,26 +80,37 @@ public:
worklet.Run(it, particles, MaxSteps);
result = ParticleAdvectionResult<ParticleType>(particles);
}
template <typename IntegratorType, typename ParticleType, typename ParticleStorage>
ParticleAdvectionResult<ParticleType> Run(
const IntegratorType& it,
vtkm::cont::ArrayHandle<ParticleType, ParticleStorage>& particles,
vtkm::Id MaxSteps)
*/
template <typename IntegratorType,
typename ParticleType,
typename ParticleStorage,
typename TerminationType,
typename AnalysisType>
void Run(const IntegratorType& it,
vtkm::cont::ArrayHandle<ParticleType, ParticleStorage>& particles,
const TerminationType& termination,
AnalysisType& analysis)
{
vtkm::worklet::flow::ParticleAdvectionWorklet<IntegratorType, ParticleType> worklet;
worklet.Run(it, particles, MaxSteps);
return ParticleAdvectionResult<ParticleType>(particles);
vtkm::worklet::flow::
ParticleAdvectionWorklet<IntegratorType, ParticleType, TerminationType, AnalysisType>
worklet;
worklet.Run(it, particles, termination, analysis);
//return ParticleAdvectionResult<ParticleType>(particles);
}
template <typename IntegratorType, typename ParticleType, typename PointStorage>
ParticleAdvectionResult<ParticleType> Run(
const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec3f, PointStorage>& points,
vtkm::Id MaxSteps)
template <typename IntegratorType,
typename ParticleType,
typename PointStorage,
typename TerminationType,
typename AnalysisType>
void Run(const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec3f, PointStorage>& points,
const TerminationType& termination,
AnalysisType& analysis)
{
vtkm::worklet::flow::ParticleAdvectionWorklet<IntegratorType, ParticleType> worklet;
vtkm::worklet::flow::
ParticleAdvectionWorklet<IntegratorType, ParticleType, TerminationType, AnalysisType>
worklet;
vtkm::cont::ArrayHandle<ParticleType> particles;
vtkm::cont::ArrayHandle<vtkm::Id> step, ids;
@ -117,11 +128,11 @@ public:
vtkm::cont::ArrayCopy(id, ids);
invoke(detail::CopyToParticle{}, points, ids, time, step, particles);
worklet.Run(it, particles, MaxSteps);
return ParticleAdvectionResult<ParticleType>(particles);
worklet.Run(it, particles, termination, analysis);
//return ParticleAdvectionResult<ParticleType>(particles);
}
};
/*
template <typename ParticleType>
struct StreamlineResult
{
@ -167,7 +178,7 @@ public:
return StreamlineResult<ParticleType>(particles, positions, polyLines);
}
};
*/
}
}
} // vtkm::worklet::flow

71
vtkm/filter/flow/worklet/ParticleAdvectionWorklets.h
View File

@ -15,6 +15,7 @@
#include <vtkm/cont/CellSetExplicit.h>
#include <vtkm/cont/ConvertNumComponentsToOffsets.h>
#include <vtkm/cont/ExecutionObjectBase.h>
#include <vtkm/cont/Invoker.h>
#include <vtkm/Particle.h>
#include <vtkm/filter/flow/worklet/Particles.h>
@ -34,18 +35,26 @@ namespace flow
class ParticleAdvectWorklet : public vtkm::worklet::WorkletMapField
{
public:
using ControlSignature = void(FieldIn idx,
ExecObject integrator,
ExecObject integralCurve,
FieldIn maxSteps);
using ExecutionSignature = void(_1 idx, _2 integrator, _3 integralCurve, _4 maxSteps);
VTKM_EXEC_CONT
ParticleAdvectWorklet()
: PushOutOfBounds(true)
{
}
VTKM_EXEC_CONT
ParticleAdvectWorklet(bool pushOutOfBounds)
: PushOutOfBounds(pushOutOfBounds)
{
}
using ControlSignature = void(FieldIn idx, ExecObject integrator, ExecObject integralCurve);
using ExecutionSignature = void(_1 idx, _2 integrator, _3 integralCurve);
using InputDomain = _1;
template <typename IntegratorType, typename IntegralCurveType>
VTKM_EXEC void operator()(const vtkm::Id& idx,
const IntegratorType& integrator,
IntegralCurveType& integralCurve,
const vtkm::Id& maxSteps) const
IntegralCurveType& integralCurve) const
{
auto particle = integralCurve.GetParticle(idx);
vtkm::FloatDefault time = particle.GetTime();
@ -69,7 +78,7 @@ public:
//We can't take a step inside spatial boundary.
//Try and take a step just past the boundary.
else if (status.CheckSpatialBounds())
else if (status.CheckSpatialBounds() && this->PushOutOfBounds)
{
status = integrator.SmallStep(particle, time, outpos);
if (status.CheckOk())
@ -78,16 +87,22 @@ public:
tookAnySteps = true;
}
}
integralCurve.StatusUpdate(idx, status, maxSteps);
integralCurve.StatusUpdate(idx, status);
} while (integralCurve.CanContinue(idx));
//Mark if any steps taken
integralCurve.UpdateTookSteps(idx, tookAnySteps);
}
private:
bool PushOutOfBounds;
};
template <typename IntegratorType, typename ParticleType>
template <typename IntegratorType,
typename ParticleType,
typename TerminationType,
typename AnalysisType>
class ParticleAdvectionWorklet
{
public:
@ -97,19 +112,17 @@ public:
void Run(const IntegratorType& integrator,
vtkm::cont::ArrayHandle<ParticleType>& particles,
vtkm::Id& MaxSteps)
const TerminationType& termination,
AnalysisType& analysis)
{
using ParticleAdvectWorkletType = vtkm::worklet::flow::ParticleAdvectWorklet;
using ParticleWorkletDispatchType =
typename vtkm::worklet::DispatcherMapField<ParticleAdvectWorkletType>;
using ParticleArrayType = vtkm::worklet::flow::Particles<ParticleType>;
using ParticleArrayType =
vtkm::worklet::flow::Particles<ParticleType, TerminationType, AnalysisType>;
vtkm::Id numSeeds = static_cast<vtkm::Id>(particles.GetNumberOfValues());
//Create and invoke the particle advection.
vtkm::cont::ArrayHandleConstant<vtkm::Id> maxSteps(MaxSteps, numSeeds);
//vtkm::cont::ArrayHandleConstant<vtkm::Id> maxSteps(MaxSteps, numSeeds);
vtkm::cont::ArrayHandleIndex idxArray(numSeeds);
// TODO: The particle advection sometimes behaves incorrectly on CUDA if the stack size
// is not changed thusly. This is concerning as the compiler should be able to determine
// staticly the required stack depth. What is even more concerning is that the runtime
@ -121,15 +134,23 @@ public:
(void)stack;
#endif // VTKM_CUDA
ParticleArrayType particlesObj(particles, MaxSteps);
// Initialize all the pre-requisites needed to start analysis
// It's based on the existing properties of the particles,
// for e.g. the number of steps they've already taken
analysis.InitializeAnalysis(particles);
//Invoke particle advection worklet
ParticleWorkletDispatchType particleWorkletDispatch;
ParticleArrayType particlesObj(particles, termination, analysis);
particleWorkletDispatch.Invoke(idxArray, integrator, particlesObj, maxSteps);
vtkm::worklet::flow::ParticleAdvectWorklet worklet(analysis.SupportPushOutOfBounds());
vtkm::cont::Invoker invoker;
invoker(worklet, idxArray, integrator, particlesObj);
// Finalize the analysis and clear intermittant arrays.
analysis.FinalizeAnalysis(particles);
}
};
/*
namespace detail
{
class GetSteps : public vtkm::worklet::WorkletMapField
@ -166,8 +187,8 @@ public:
}
};
} // namespace detail
*/
/*
template <typename IntegratorType, typename ParticleType>
class StreamlineWorklet
{
@ -231,7 +252,7 @@ public:
polyLines.Fill(positions.GetNumberOfValues(), cellTypes, connectivity, offsets);
}
};
*/
}
}
} // namespace vtkm::worklet::flow

214
vtkm/filter/flow/worklet/Particles.h
View File

@ -14,6 +14,9 @@
#include <vtkm/cont/ExecutionObjectBase.h>
#include <vtkm/filter/flow/worklet/IntegratorStatus.h>
#include <vtkm/filter/flow/worklet/Analysis.h>
#include <vtkm/filter/flow/worklet/Termination.h>
namespace vtkm
{
namespace worklet
@ -21,24 +24,27 @@ namespace worklet
namespace flow
{
template <typename ParticleType>
template <typename ParticleType, typename TerminationType, typename AnalysisType>
class ParticleExecutionObject
{
public:
VTKM_EXEC_CONT
ParticleExecutionObject()
: Particles()
, MaxSteps(0)
, Termination()
, Analysis()
{
}
ParticleExecutionObject(vtkm::cont::ArrayHandle<ParticleType> particleArray,
vtkm::Id maxSteps,
const TerminationType& termination,
const AnalysisType& analysis,
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token)
: Termination(termination)
, Analysis(analysis)
{
Particles = particleArray.PrepareForInPlace(device, token);
MaxSteps = maxSteps;
}
VTKM_EXEC
@ -57,19 +63,15 @@ public:
newParticle.SetPosition(pt);
newParticle.SetTime(time);
newParticle.SetNumberOfSteps(particle.GetNumberOfSteps() + 1);
this->Analysis.Analyze(idx, particle, newParticle);
this->Particles.Set(idx, newParticle);
}
VTKM_EXEC
void StatusUpdate(const vtkm::Id& idx,
const vtkm::worklet::flow::IntegratorStatus& status,
vtkm::Id maxSteps)
void StatusUpdate(const vtkm::Id& idx, const vtkm::worklet::flow::IntegratorStatus& status)
{
ParticleType p(this->GetParticle(idx));
if (p.GetNumberOfSteps() == maxSteps)
p.GetStatus().SetTerminate();
if (status.CheckFail())
p.GetStatus().SetFail();
if (status.CheckSpatialBounds())
@ -88,7 +90,13 @@ public:
}
VTKM_EXEC
bool CanContinue(const vtkm::Id& idx) { return this->GetParticle(idx).GetStatus().CanContinue(); }
bool CanContinue(const vtkm::Id& idx)
{
ParticleType particle(this->GetParticle(idx));
auto terminate = this->Termination.CheckTermination(particle);
this->Particles.Set(idx, particle);
return terminate;
}
VTKM_EXEC
void UpdateTookSteps(const vtkm::Id& idx, bool val)
@ -103,176 +111,48 @@ public:
protected:
using ParticlePortal = typename vtkm::cont::ArrayHandle<ParticleType>::WritePortalType;
ParticlePortal Particles;
vtkm::Id MaxSteps;
TerminationType Termination;
AnalysisType Analysis;
};
template <typename ParticleType>
template <typename ParticleType, typename TerminationType, typename AnalysisType>
class Particles : public vtkm::cont::ExecutionObjectBase
{
protected:
vtkm::cont::ArrayHandle<ParticleType> ParticleArray;
TerminationType Termination;
AnalysisType Analysis;
public:
VTKM_CONT vtkm::worklet::flow::ParticleExecutionObject<ParticleType> PrepareForExecution(
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
{
return vtkm::worklet::flow::ParticleExecutionObject<ParticleType>(
this->ParticleArray, this->MaxSteps, device, token);
}
VTKM_CONT
Particles() = default;
VTKM_CONT
Particles(vtkm::cont::ArrayHandle<ParticleType>& pArray, vtkm::Id& maxSteps)
Particles(vtkm::cont::ArrayHandle<ParticleType>& pArray,
TerminationType termination,
AnalysisType analysis)
: ParticleArray(pArray)
, MaxSteps(maxSteps)
, Termination(termination)
, Analysis(analysis)
{
}
Particles() {}
protected:
vtkm::cont::ArrayHandle<ParticleType> ParticleArray;
vtkm::Id MaxSteps;
VTKM_CONT auto PrepareForExecution(vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
-> vtkm::worklet::flow::ParticleExecutionObject<
ParticleType,
decltype(this->Termination.PrepareForExecution(device, token)),
decltype(this->Analysis.PrepareForExecution(device, token))>
{
auto termination = this->Termination.PrepareForExecution(device, token);
auto analysis = this->Analysis.PrepareForExecution(device, token);
return vtkm::worklet::flow::
ParticleExecutionObject<ParticleType, decltype(termination), decltype(analysis)>(
this->ParticleArray, termination, analysis, device, token);
}
};
template <typename ParticleType>
class StateRecordingParticleExecutionObject : public ParticleExecutionObject<ParticleType>
{
public:
VTKM_EXEC_CONT
StateRecordingParticleExecutionObject()
: ParticleExecutionObject<ParticleType>()
, History()
, Length(0)
, StepCount()
, ValidPoint()
{
}
StateRecordingParticleExecutionObject(vtkm::cont::ArrayHandle<ParticleType> pArray,
vtkm::cont::ArrayHandle<vtkm::Vec3f> historyArray,
vtkm::cont::ArrayHandle<vtkm::Id> validPointArray,
vtkm::cont::ArrayHandle<vtkm::Id> stepCountArray,
vtkm::Id maxSteps,
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token)
: ParticleExecutionObject<ParticleType>(pArray, maxSteps, device, token)
, Length(maxSteps + 1)
{
vtkm::Id numPos = pArray.GetNumberOfValues();
History = historyArray.PrepareForOutput(numPos * Length, device, token);
ValidPoint = validPointArray.PrepareForInPlace(device, token);
StepCount = stepCountArray.PrepareForInPlace(device, token);
}
VTKM_EXEC
void PreStepUpdate(const vtkm::Id& idx)
{
ParticleType p = this->ParticleExecutionObject<ParticleType>::GetParticle(idx);
if (this->StepCount.Get(idx) == 0)
{
vtkm::Id loc = idx * Length;
this->History.Set(loc, p.GetPosition());
this->ValidPoint.Set(loc, 1);
this->StepCount.Set(idx, 1);
}
}
VTKM_EXEC
void StepUpdate(const vtkm::Id& idx,
const ParticleType& particle,
vtkm::FloatDefault time,
const vtkm::Vec3f& pt)
{
this->ParticleExecutionObject<ParticleType>::StepUpdate(idx, particle, time, pt);
//local step count.
vtkm::Id stepCount = this->StepCount.Get(idx);
vtkm::Id loc = idx * Length + stepCount;
this->History.Set(loc, pt);
this->ValidPoint.Set(loc, 1);
this->StepCount.Set(idx, stepCount + 1);
}
protected:
using IdPortal = typename vtkm::cont::ArrayHandle<vtkm::Id>::WritePortalType;
using HistoryPortal = typename vtkm::cont::ArrayHandle<vtkm::Vec3f>::WritePortalType;
HistoryPortal History;
vtkm::Id Length;
IdPortal StepCount;
IdPortal ValidPoint;
};
template <typename ParticleType>
class StateRecordingParticles : vtkm::cont::ExecutionObjectBase
{
public:
//Helper functor for compacting history
struct IsOne
{
template <typename T>
VTKM_EXEC_CONT bool operator()(const T& x) const
{
return x == T(1);
}
};
VTKM_CONT vtkm::worklet::flow::StateRecordingParticleExecutionObject<ParticleType>
PrepareForExecution(vtkm::cont::DeviceAdapterId device, vtkm::cont::Token& token) const
{
return vtkm::worklet::flow::StateRecordingParticleExecutionObject<ParticleType>(
this->ParticleArray,
this->HistoryArray,
this->ValidPointArray,
this->StepCountArray,
this->MaxSteps,
device,
token);
}
VTKM_CONT
StateRecordingParticles(vtkm::cont::ArrayHandle<ParticleType>& pArray, const vtkm::Id& maxSteps)
: MaxSteps(maxSteps)
, ParticleArray(pArray)
{
vtkm::Id numParticles = static_cast<vtkm::Id>(pArray.GetNumberOfValues());
//Create ValidPointArray initialized to zero.
vtkm::cont::ArrayHandleConstant<vtkm::Id> tmp(0, (this->MaxSteps + 1) * numParticles);
vtkm::cont::ArrayCopy(tmp, this->ValidPointArray);
//Create StepCountArray initialized to zero.
vtkm::cont::ArrayHandleConstant<vtkm::Id> tmp2(0, numParticles);
vtkm::cont::ArrayCopy(tmp2, this->StepCountArray);
}
VTKM_CONT
StateRecordingParticles(vtkm::cont::ArrayHandle<ParticleType>& pArray,
vtkm::cont::ArrayHandle<vtkm::Vec3f>& historyArray,
vtkm::cont::ArrayHandle<vtkm::Id>& validPointArray,
vtkm::Id& maxSteps)
{
ParticleArray = pArray;
HistoryArray = historyArray;
ValidPointArray = validPointArray;
MaxSteps = maxSteps;
}
VTKM_CONT
void GetCompactedHistory(vtkm::cont::ArrayHandle<vtkm::Vec3f>& positions)
{
vtkm::cont::Algorithm::CopyIf(this->HistoryArray, this->ValidPointArray, positions, IsOne());
}
protected:
vtkm::cont::ArrayHandle<vtkm::Vec3f> HistoryArray;
vtkm::Id MaxSteps;
vtkm::cont::ArrayHandle<ParticleType> ParticleArray;
vtkm::cont::ArrayHandle<vtkm::Id> StepCountArray;
vtkm::cont::ArrayHandle<vtkm::Id> ValidPointArray;
};
}
}
} //vtkm::worklet::flow

2
vtkm/filter/flow/worklet/RK4Integrator.h
View File

@ -34,7 +34,7 @@ public:
}
template <typename Particle>
VTKM_EXEC IntegratorStatus CheckStep(Particle& particle,
VTKM_EXEC IntegratorStatus CheckStep(const Particle& particle,
vtkm::FloatDefault stepLength,
vtkm::Vec3f& velocity) const
{

173
vtkm/filter/flow/worklet/Termination.h Normal file
View File

@ -0,0 +1,173 @@
//=============================================================================
//
// 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 vtkm_worklet_particleadvection_termination
#define vtkm_worklet_particleadvection_termination
#include <vtkm/Types.h>
#include <vtkm/cont/ExecutionObjectBase.h>
namespace vtkm
{
namespace worklet
{
namespace flow
{
class NormalTerminationExec
{
public:
VTKM_EXEC_CONT
NormalTerminationExec()
: MaxSteps(0)
{
}
VTKM_EXEC_CONT
NormalTerminationExec(vtkm::Id maxSteps)
: MaxSteps(maxSteps)
{
}
template <typename ParticleType>
VTKM_EXEC bool CheckTermination(ParticleType& particle) const
{
/// Checks particle properties to make a decision for termination
/// -- Check if the particle is out of spatial boundaries
/// -- Check if the particle has reached the maximum number of steps
/// -- Check if the particle is in a zero velocity region
auto& status = particle.GetStatus();
if (particle.GetNumberOfSteps() == this->MaxSteps)
{
status.SetTerminate();
particle.SetStatus(status);
}
bool terminate = status.CheckOk() && !status.CheckTerminate() && !status.CheckSpatialBounds() &&
!status.CheckTemporalBounds() && !status.CheckInGhostCell() && !status.CheckZeroVelocity();
return terminate;
}
private:
vtkm::Id MaxSteps;
};
class NormalTermination : public vtkm::cont::ExecutionObjectBase
{
public:
VTKM_CONT
NormalTermination()
: MaxSteps(0)
{
}
VTKM_CONT
NormalTermination(vtkm::Id maxSteps)
: MaxSteps(maxSteps)
{
}
VTKM_CONT
vtkm::Id AllocationSize() const { return this->MaxSteps; }
VTKM_CONT vtkm::worklet::flow::NormalTerminationExec PrepareForExecution(
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
{
(void)device;
(void)token;
return vtkm::worklet::flow::NormalTerminationExec(MaxSteps);
}
private:
vtkm::Id MaxSteps;
};
class PoincareTerminationExec
{
public:
VTKM_EXEC_CONT
PoincareTerminationExec()
: MaxSteps(0)
, MaxPunctures(0)
{
}
VTKM_EXEC_CONT
PoincareTerminationExec(vtkm::Id maxSteps, vtkm::Id maxPunctures)
: MaxSteps(maxSteps)
, MaxPunctures(maxPunctures)
{
}
template <typename ParticleType>
VTKM_EXEC bool CheckTermination(ParticleType& particle) const
{
/// Checks particle properties to make a decision for termination
/// -- Check if the particle is out of spatial boundaries
/// -- Check if the particle has reached the maximum number of steps
/// -- Check if the particle is in a zero velocity region
auto& status = particle.GetStatus();
if (particle.GetNumberOfSteps() >= this->MaxSteps ||
particle.GetNumberOfPunctures() >= this->MaxPunctures)
{
status.SetTerminate();
particle.SetStatus(status);
}
bool terminate = status.CheckOk() && !status.CheckTerminate() && !status.CheckSpatialBounds() &&
!status.CheckTemporalBounds() && !status.CheckInGhostCell() && !status.CheckZeroVelocity();
return terminate;
}
private:
vtkm::Id MaxSteps;
vtkm::Id MaxPunctures;
};
class PoincareTermination : public vtkm::cont::ExecutionObjectBase
{
public:
VTKM_CONT
PoincareTermination()
: MaxSteps(0)
, MaxPunctures(0)
{
}
VTKM_CONT
PoincareTermination(vtkm::Id maxSteps, vtkm::Id maxPunctures)
: MaxSteps(maxSteps)
, MaxPunctures(maxPunctures)
{
}
VTKM_CONT
vtkm::Id AllocationSize() const { return this->MaxPunctures; }
VTKM_CONT vtkm::worklet::flow::PoincareTerminationExec PrepareForExecution(
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
{
(void)device;
(void)token;
return vtkm::worklet::flow::PoincareTerminationExec(MaxSteps, MaxPunctures);
}
private:
vtkm::Id MaxSteps;
vtkm::Id MaxPunctures;
};
} // namespace particleadvection
} // namespace worklet
} // namespace vtkm
#endif

400
vtkm/filter/flow/worklet/XGCField.h Normal file
View File

@ -0,0 +1,400 @@
//============================================================================
// 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 vtkm_filter_flow_worklet_XGCField_h
#define vtkm_filter_flow_worklet_XGCField_h
#include <vtkm/Geometry.h>
#include <vtkm/Matrix.h>
#include <vtkm/Types.h>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/VecVariable.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ExecutionObjectBase.h>
#include <vtkm/exec/CellInterpolate.h>
#include <vtkm/filter/flow/worklet/XGCHelper.h>
namespace vtkm
{
namespace worklet
{
namespace flow
{
template <typename FieldComponentType, typename FieldVecType>
class ExecutionXGCField
{
public:
using ComponentPortalType = typename FieldComponentType::ReadPortalType;
using VecPortalType = typename FieldVecType::ReadPortalType;
using Association = vtkm::cont::Field::Association;
using DelegateToField = std::true_type;
using Ray3f = vtkm::Ray<vtkm::FloatDefault, 3, true>;
VTKM_CONT
ExecutionXGCField(const FieldComponentType& as_ff,
const FieldVecType& _dAs_ff_rzp,
const FieldComponentType& coeff_1D,
const FieldComponentType& coeff_2D,
const FieldVecType& _B_rzp,
const FieldComponentType& psi,
const XGCParams& params,
vtkm::FloatDefault period,
bool useDeltaBScale,
vtkm::FloatDefault deltaBScale,
bool useBScale,
vtkm::FloatDefault bScale,
const Association assoc,
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token)
: As_ff(as_ff.PrepareForInput(device, token))
, dAs_ff_rzp(_dAs_ff_rzp.PrepareForInput(device, token))
, Coeff_1D(coeff_1D.PrepareForInput(device, token))
, Coeff_2D(coeff_2D.PrepareForInput(device, token))
, Psi(psi.PrepareForInput(device, token))
, B_rzp(_B_rzp.PrepareForInput(device, token))
, Params(params)
, Period(period)
, UseDeltaBScale(useDeltaBScale)
, DeltaBScale(deltaBScale)
, UseBScale(useBScale)
, BScale(bScale)
, Assoc(assoc)
{
}
private:
template <typename LocatorType, typename InterpolationHelper>
VTKM_EXEC bool PtLoc(const vtkm::Vec3f& ptRZ,
const LocatorType& locator,
const InterpolationHelper& helper,
vtkm::Vec3f& param,
vtkm::Vec<vtkm::Id, 3>& vIds) const
{
vtkm::Id cellId;
vtkm::ErrorCode status;
status = locator.FindCell(ptRZ, cellId, param);
if (status != vtkm::ErrorCode::Success)
{
return false;
}
vtkm::UInt8 cellShape;
vtkm::IdComponent nVerts;
vtkm::VecVariable<vtkm::Id, 8> ptIndices;
helper.GetCellInfo(cellId, cellShape, nVerts, ptIndices);
vIds[0] = ptIndices[0];
vIds[1] = ptIndices[1];
vIds[2] = ptIndices[2];
return true;
}
public:
VTKM_EXEC Association GetAssociation() const { return this->Assoc; }
VTKM_EXEC void GetValue(const vtkm::Id vtkmNotUsed(cellId),
vtkm::VecVariable<vtkm::Vec3f, 2>& vtkmNotUsed(value)) const
{
//TODO Raise Error : XGC Field should not allow this path
}
VTKM_EXEC void GetValue(const vtkm::VecVariable<vtkm::Id, 8>& vtkmNotUsed(indices),
const vtkm::Id vtkmNotUsed(vertices),
const vtkm::Vec3f& vtkmNotUsed(parametric),
const vtkm::UInt8 vtkmNotUsed(cellShape),
vtkm::VecVariable<vtkm::Vec3f, 2>& vtkmNotUsed(value)) const
{
//TODO Raise Error : XGC Field should not allow this path
}
template <typename Point, typename Locator, typename Helper>
VTKM_EXEC bool GetValue(const Point& point,
const vtkm::FloatDefault& time,
vtkm::VecVariable<Point, 2>& out,
const Locator& locator,
const Helper& helper) const
{
// not used
(void)time;
auto R = point[0];
auto Phi = point[1];
auto Z = point[2];
ParticleMetaData metadata;
if (!HighOrderB(point, metadata, this->Params, this->Coeff_1D, this->Coeff_2D))
return false;
vtkm::Id planeIdx0, planeIdx1, numRevs;
vtkm::FloatDefault phiN, Phi0, Phi1, T;
GetPlaneIdx(
Phi, this->Period, this->Params, phiN, planeIdx0, planeIdx1, Phi0, Phi1, numRevs, T);
vtkm::Vec3f B0_rpz(metadata.B0_rzp[0], metadata.B0_rzp[2], metadata.B0_rzp[1]);
vtkm::Vec3f ff_pt_rpz;
CalcFieldFollowingPt(
{ R, phiN, Z }, B0_rpz, Phi0, Phi1, this->Params, this->Coeff_1D, this->Coeff_2D, ff_pt_rpz);
//Now, interpolate between Phi_i and Phi_i+1
vtkm::FloatDefault T01 = (phiN - Phi0) / (Phi1 - Phi0);
vtkm::FloatDefault T10 = 1.0f - T01;
//Get vec at Phi0 and Phi1.
//x_ff is in rzp
//vtkm::Vec3f x_ff_rzp(ptOnMidPlane_rpz[0], ptOnMidPlane_rpz[2], 0);
vtkm::Vec3f x_ff_rzp(ff_pt_rpz[0], ff_pt_rpz[2], 0);
// Offsets into As and DAs to jump to the right plane.
int offsets[2];
offsets[0] = planeIdx0 * this->Params.NumNodes * 2;
offsets[1] = planeIdx0 * this->Params.NumNodes * 2 + this->Params.NumNodes;
const vtkm::FloatDefault basis = 0.0f;
//auto B0_R = B0_rzp[0];
//auto B0_Z = B0_rzp[1];
//auto x_ff_R = x_ff_rzp[0];
//auto x_ff_Z = x_ff_rzp[1];
//gradPsi: pt on mid plane? (question)
//dPsi/dR = B0_Z * R
//dPsi/dZ = -B0_R * R;
//vtkm::Vec3f gradPsi_rzp(B0_Z * x_ff_R, -B0_R * x_ff_R, 0);
//use high order...
vtkm::Vec3f gradPsi_rzp = metadata.gradPsi_rzp;
vtkm::FloatDefault gammaPsi = 1.0f / vtkm::Magnitude(gradPsi_rzp);
vtkm::Vec2f rvec(0, 0), zvec(0, 0);
rvec[0] = basis + (1.0 - basis) * gammaPsi * gradPsi_rzp[0];
rvec[1] = (1.0 - basis) * gammaPsi * gradPsi_rzp[1];
zvec[0] = (1.0 - basis) * gammaPsi * (-gradPsi_rzp[1]);
zvec[1] = basis + (1.0 - basis) * gammaPsi * gradPsi_rzp[0];
//Get the vectors in the ff coordinates.
//auto dAs_ff_rzp = EvalVector(ds, locator, {x_ff_rzp, x_ff_rzp}, "dAs_ff_rzp", offsets);
//auto dAs_ff0_rzp = dAs_ff_rzp[0];
//auto dAs_ff1_rzp = dAs_ff_rzp[1];
vtkm::Vec3f x_ff_param;
vtkm::Vec<vtkm::Id, 3> x_ff_vids;
//Hotspot..
if (!this->PtLoc(x_ff_rzp, locator, helper, x_ff_param, x_ff_vids))
return false;
// Eval actual values for DAsPhiFF from the triangular planes
auto dAs_ff0_rzp = Eval(this->dAs_ff_rzp, offsets[0], x_ff_param, x_ff_vids);
auto dAs_ff1_rzp = Eval(this->dAs_ff_rzp, offsets[1], x_ff_param, x_ff_vids);
vtkm::FloatDefault wphi[2] = { T10, T01 }; //{T01, T10};
vtkm::Vec3f gradAs_rpz;
//vec.r = wphi[0]*( rvec[0]*V.r[0] + zvec[0]*V.z[0]) +
// wphi[1]*( rvec[0]*V.r[1] + zvec[0]*v.z[1]);
//vec.p = wphi[0]*V.phi[0] +
// whpi[1]*V.phi[1];
//vec.z = wphi[0]*( rvec[1]*V.r[0] + zvec[1]*V.z[0]) +
// wphi[1]*( rvec[1]*V.r[1] + zvec[1]*V.Z[1]);
gradAs_rpz[0] = wphi[0] * (rvec[0] * dAs_ff0_rzp[0] + zvec[0] * dAs_ff0_rzp[1]) +
wphi[1] * (rvec[0] * dAs_ff1_rzp[0] + zvec[0] * dAs_ff1_rzp[1]);
gradAs_rpz[1] = wphi[0] * dAs_ff0_rzp[2] + wphi[1] * dAs_ff1_rzp[2];
gradAs_rpz[2] = wphi[0] * (rvec[1] * dAs_ff0_rzp[0] + zvec[1] * dAs_ff0_rzp[1]) +
wphi[1] * (rvec[1] * dAs_ff1_rzp[0] + zvec[1] * dAs_ff1_rzp[1]);
vtkm::FloatDefault BMag = vtkm::Magnitude(metadata.B0_rzp);
//project using bfield.
//gradAs.Phi = (gradAs.Phi * BMag - gradAs.R*B0_pos.R - gradAs.Z*B0_pos.Z) / B0_pos.Phi
gradAs_rpz[1] = (gradAs_rpz[1] * BMag - gradAs_rpz[0] * metadata.B0_rzp[0] -
gradAs_rpz[2] * metadata.B0_rzp[1]) /
metadata.B0_rzp[2];
//deltaB = AsCurl(bhat) + gradAs x bhat.
//std::vector<int> off = {planeIdx0*this->NumNodes};
//vtkm::Vec3f AsCurl_bhat_rzp = EvalVector(ds, locator, {x_ff_rzp}, "AsCurlBHat_RZP", off)[0];
//auto AsCurl_bhat_rzp = this->EvalV(AsCurlBHat_RZP, 0, x_ff_vids, x_ff_param);
//
auto As_ff0 = Eval(this->As_ff, offsets[0], x_ff_param, x_ff_vids);
auto As_ff1 = Eval(this->As_ff, offsets[1], x_ff_param, x_ff_vids);
// Interpolated value of As on plane 10 & 1.
// wPhi is distance of the plane from the point.
vtkm::FloatDefault As = wphi[0] * As_ff0 + wphi[1] * As_ff1;
auto AsCurl_bhat_rzp = As * metadata.curl_nb_rzp;
auto bhat_rzp = vtkm::Normal(metadata.B0_rzp);
vtkm::Vec3f gradAs_rzp(gradAs_rpz[0], gradAs_rpz[2], gradAs_rpz[1]);
vtkm::Vec3f deltaB_rzp = AsCurl_bhat_rzp + vtkm::Cross(gradAs_rzp, bhat_rzp);
if (this->UseDeltaBScale)
deltaB_rzp = deltaB_rzp * this->DeltaBScale;
vtkm::Vec3f B0_rzp = metadata.B0_rzp;
if (this->UseBScale)
B0_rzp = B0_rzp * this->BScale;
deltaB_rzp[2] /= R;
B0_rzp[2] /= R;
vtkm::Vec3f vec_rzp = B0_rzp + deltaB_rzp;
vtkm::Vec3f vec_rpz(vec_rzp[0], vec_rzp[2], vec_rzp[1]);
out = vtkm::make_Vec(vec_rpz);
return true;
}
private:
ComponentPortalType As_ff;
VecPortalType dAs_ff_rzp;
ComponentPortalType Coeff_1D;
ComponentPortalType Coeff_2D;
ComponentPortalType Psi;
VecPortalType B_rzp;
vtkm::FloatDefault Period;
bool UseDeltaBScale;
vtkm::FloatDefault DeltaBScale = 1.0;
bool UseBScale;
vtkm::FloatDefault BScale = 1.0;
Association Assoc;
XGCParams Params;
};
template <typename FieldComponentType, typename FieldVecType>
class XGCField : public vtkm::cont::ExecutionObjectBase
{
public:
using ExecutionType = ExecutionXGCField<FieldComponentType, FieldVecType>;
using Association = vtkm::cont::Field::Association;
VTKM_CONT
XGCField() = default;
VTKM_CONT
XGCField(const FieldComponentType& as_ff,
const FieldVecType& _dAs_ff_rzp,
const FieldComponentType& coeff_1D,
const FieldComponentType& coeff_2D,
const FieldVecType& _B_rzp,
const FieldComponentType& psi,
const XGCParams& params,
vtkm::FloatDefault period,
bool useDeltaBScale,
vtkm::FloatDefault deltaBScale,
bool useBScale,
vtkm::FloatDefault bScale)
: As_ff(as_ff)
, dAs_ff_rzp(_dAs_ff_rzp)
, Coeff_1D(coeff_1D)
, Coeff_2D(coeff_2D)
, Psi(psi)
, B_rzp(_B_rzp)
, Params(params)
, Period(period)
, UseDeltaBScale(useDeltaBScale)
, DeltaBScale(deltaBScale)
, UseBScale(useBScale)
, BScale(bScale)
, Assoc(vtkm::cont::Field::Association::Points)
{
}
VTKM_CONT
XGCField(const FieldComponentType& as_ff,
const FieldVecType& _dAs_ff_rzp,
const FieldComponentType& coeff_1D,
const FieldComponentType& coeff_2D,
const FieldVecType& _B_rzp,
const FieldComponentType& psi,
const XGCParams& params,
vtkm::FloatDefault period,
bool useDeltaBScale,
vtkm::FloatDefault deltaBScale,
bool useBScale,
vtkm::FloatDefault bScale,
const Association assoc)
: As_ff(as_ff)
, dAs_ff_rzp(_dAs_ff_rzp)
, Coeff_1D(coeff_1D)
, Coeff_2D(coeff_2D)
, Psi(psi)
, B_rzp(_B_rzp)
, Params(params)
, Period(period)
, UseDeltaBScale(useDeltaBScale)
, DeltaBScale(deltaBScale)
, UseBScale(useBScale)
, BScale(bScale)
, Assoc(assoc)
{
if (assoc != Association::Points && assoc != Association::Cells)
throw("Unsupported field association");
}
VTKM_CONT
const ExecutionType PrepareForExecution(vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
{
return ExecutionType(this->As_ff,
this->dAs_ff_rzp,
this->Coeff_1D,
this->Coeff_2D,
this->B_rzp,
this->Psi,
this->Params,
this->Period,
this->UseDeltaBScale,
this->DeltaBScale,
this->UseBScale,
this->BScale,
this->Assoc,
device,
token);
}
private:
// For turbulance in the magnetic field (B is constant, but this emulates deltaB)
FieldComponentType As_ff;
FieldVecType dAs_ff_rzp;
// For Higher order interpolation
// on RZ uniform grid, not the triangle grid.
FieldComponentType Coeff_1D;
FieldComponentType Coeff_2D;
// Constant magnetic field on the triangle mesh
FieldComponentType Psi;
// B field at each vertex of the triangular field
FieldVecType B_rzp;
XGCParams Params;
vtkm::FloatDefault Period;
bool UseDeltaBScale;
vtkm::FloatDefault DeltaBScale = 1.0;
bool UseBScale;
vtkm::FloatDefault BScale = 1.0;
Association Assoc;
};
}
}
} //vtkm::worklet::flow
#endif //vtkm_filter_flow_worklet_XGCField_h

507
vtkm/filter/flow/worklet/XGCHelper.h Normal file
View File

@ -0,0 +1,507 @@
//============================================================================
// 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 vtkm_filter_flow_worklet_XGCHelper_h
#define vtkm_filter_flow_worklet_XGCHelper_h
#include <vtkm/Geometry.h>
#include <vtkm/Matrix.h>
namespace vtkm
{
namespace worklet
{
namespace flow
{
struct XGCParams
{
vtkm::Id NumPlanes = -1;
vtkm::Id NumNodes = -1;
vtkm::Id NumTri = -1;
vtkm::FloatDefault Period;
vtkm::FloatDefault dPhi;
int ncoeff_r, ncoeff_z, ncoeff_psi;
vtkm::FloatDefault min_psi, max_psi;
vtkm::FloatDefault sml_bp_sign = -1.0f;
vtkm::FloatDefault eq_axis_r = -1, eq_axis_z = -1, eq_x_psi = -1, eq_x_r = -1, eq_x_z = -1;
vtkm::FloatDefault eq_min_r = -1, eq_max_r = -1;
vtkm::FloatDefault eq_min_z = -1, eq_max_z = 1;
vtkm::FloatDefault psi_min = -1, psi_max = -1;
vtkm::Id eq_mr = -1, eq_mz = -1;
vtkm::Id eq_mpsi = -1;
vtkm::Id sml_wedge_n = -1;
vtkm::FloatDefault itp_min_psi = -1, itp_max_psi = -1;
vtkm::FloatDefault one_d_cub_dpsi_inv;
int nr, nz;
vtkm::FloatDefault rmin, zmin, rmax, zmax;
vtkm::FloatDefault dr, dz, dr_inv, dz_inv;
vtkm::FloatDefault eq_x2_z, eq_x2_r, eq_x2_psi;
vtkm::FloatDefault eq_x_slope, eq_x2_slope;
};
struct ParticleMetaData
{
ParticleMetaData() = default;
ParticleMetaData(const ParticleMetaData&) = default;
ParticleMetaData& operator=(const ParticleMetaData&) = default;
vtkm::Id3 PrevCell = { -1, -1, -1 };
vtkm::FloatDefault Psi = 0;
vtkm::FloatDefault dpsi_dr = 0, dpsi_dz = 0, d2psi_drdz = 0, d2psi_d2r = 0, d2psi_d2z = 0;
vtkm::Vec3f gradPsi_rzp, B0_rzp, curlB_rzp, curl_nb_rzp;
};
VTKM_EXEC
inline vtkm::Id GetIndex(const vtkm::FloatDefault& x,
const vtkm::Id& nx,
const vtkm::FloatDefault& xmin,
const vtkm::FloatDefault& dx_inv)
{
vtkm::Id idx = vtkm::Max(vtkm::Id(0), vtkm::Min(nx - 1, vtkm::Id((x - xmin) * dx_inv)));
return idx;
}
template <typename ScalarPortal>
VTKM_EXEC void EvalBicub2(const vtkm::FloatDefault& x,
const vtkm::FloatDefault& y,
const vtkm::FloatDefault& xc,
const vtkm::FloatDefault& yc,
const vtkm::Id& offset,
const ScalarPortal& Coeff_2D,
vtkm::FloatDefault& f00,
vtkm::FloatDefault& f10,
vtkm::FloatDefault& f01,
vtkm::FloatDefault& f11,
vtkm::FloatDefault& f20,
vtkm::FloatDefault& f02)
{
vtkm::FloatDefault dx = x - xc;
vtkm::FloatDefault dy = y - yc;
//fortran code.
f00 = f01 = f10 = f11 = f20 = f02 = 0.0f;
vtkm::FloatDefault xv[4] = { 1, dx, dx * dx, dx * dx * dx };
vtkm::FloatDefault yv[4] = { 1, dy, dy * dy, dy * dy * dy };
vtkm::FloatDefault fx[4] = { 0, 0, 0, 0 };
vtkm::FloatDefault dfx[4] = { 0, 0, 0, 0 };
vtkm::FloatDefault dfy[4] = { 0, 0, 0, 0 };
vtkm::FloatDefault dfx2[4] = { 0, 0, 0, 0 };
vtkm::FloatDefault dfy2[4] = { 0, 0, 0, 0 };
for (int j = 0; j < 4; j++)
{
for (int i = 0; i < 4; i++)
fx[j] = fx[j] + xv[i] * Coeff_2D.Get(offset + j * 4 + i); //acoeff[i][j];
for (int i = 1; i < 4; i++)
dfx[j] = dfx[j] +
vtkm::FloatDefault(i) * xv[i - 1] * Coeff_2D.Get(offset + j * 4 + i); //acoeff[i][j];
for (int i = 2; i < 4; i++)
dfx2[j] = dfx2[j] +
vtkm::FloatDefault(i * (i - 1)) * xv[i - 2] *
Coeff_2D.Get(offset + j * 4 + i); //acoeff[i][j];
}
for (int j = 0; j < 4; j++)
{
f00 = f00 + fx[j] * yv[j];
f10 = f10 + dfx[j] * yv[j];
f20 = f20 + dfx2[j] * yv[j];
}
for (int j = 1; j < 4; j++)
{
dfy[j] = vtkm::FloatDefault(j) * yv[j - 1];
f01 = f01 + fx[j] * dfy[j];
f11 = f11 + dfx[j] * dfy[j];
}
for (int j = 2; j < 4; j++)
{
dfy2[j] = vtkm::FloatDefault(j * (j - 1)) * yv[j - 2];
f02 = f02 + fx[j] * dfy2[j];
}
}
template <typename ScalarPortal>
VTKM_EXEC vtkm::FloatDefault I_interpol(const vtkm::FloatDefault& in_psi,
const vtkm::Id& ideriv,
const vtkm::Id& region,
const XGCParams& Params,
const ScalarPortal& Coeff_1D)
{
vtkm::FloatDefault psi;
if (region == 3)
{
psi = vtkm::Min(Params.eq_x_psi, Params.itp_max_psi);
if (ideriv != 0)
return 0.0;
}
else
{
psi = vtkm::Min(in_psi, Params.itp_max_psi);
psi = vtkm::Max(psi, Params.itp_min_psi);
}
vtkm::FloatDefault pn = psi * Params.one_d_cub_dpsi_inv;
vtkm::Id ip = floor(pn);
ip = vtkm::Min(vtkm::Max(ip, vtkm::Id(0)), vtkm::Id(Params.ncoeff_psi - 1));
vtkm::FloatDefault wp = pn - (vtkm::FloatDefault)(ip);
int idx = ip * 4;
//vtkm::FloatDefault acoef[4];
//acoef[0] = one_d_cub_acoef(ip).coeff[0];
//acoef[1] = one_d_cub_acoef(ip).coeff[1];
//acoef[2] = one_d_cub_acoef(ip).coeff[2];
//acoef[3] = one_d_cub_acoef(ip).coeff[3];
const vtkm::FloatDefault acoef[4] = {
Coeff_1D.Get(idx + 0), Coeff_1D.Get(idx + 1), Coeff_1D.Get(idx + 2), Coeff_1D.Get(idx + 3)
};
vtkm::FloatDefault iVal = 0.0;
if (ideriv == 0)
iVal = acoef[0] + (acoef[1] + (acoef[2] + acoef[3] * wp) * wp) * wp;
else if (ideriv == 1)
iVal = (acoef[1] + (2.0 * acoef[2] + 3.0 * acoef[3] * wp) * wp) * Params.one_d_cub_dpsi_inv;
return iVal * Params.sml_bp_sign;
}
VTKM_EXEC
inline bool IsRegion12(const vtkm::FloatDefault& R,
const vtkm::FloatDefault& Z,
const vtkm::FloatDefault& psi,
const XGCParams& Params)
{
constexpr vtkm::FloatDefault epsil_psi = 1e-5;
if ((psi <= (Params.eq_x_psi - epsil_psi) &&
-(R - Params.eq_x_r) * Params.eq_x_slope + (Z - Params.eq_x_z) > 0 &&
-(R - Params.eq_x2_r) * Params.eq_x2_slope + (Z - Params.eq_x2_z) < 0) ||
(psi > Params.eq_x_psi - epsil_psi && psi <= Params.eq_x2_psi - epsil_psi &&
-(R - Params.eq_x2_r) * Params.eq_x2_slope + (Z - Params.eq_x2_z) < 0) ||
psi > Params.eq_x2_psi - epsil_psi)
{
return true;
}
return false;
}
template <typename ScalarPortal>
VTKM_EXEC bool HighOrderB(const vtkm::Vec3f& pointRPZ,
ParticleMetaData& metadata,
const XGCParams& Params,
const ScalarPortal& Coeff_1D,
const ScalarPortal& Coeff_2D)
{
vtkm::FloatDefault R = pointRPZ[0], Z = pointRPZ[2];
vtkm::Id r_i = GetIndex(R, Params.nr, Params.eq_min_r, Params.dr_inv);
vtkm::Id z_i = GetIndex(Z, Params.nz, Params.zmin, Params.dz_inv);
vtkm::FloatDefault Rc = Params.eq_min_r + (vtkm::FloatDefault)(r_i)*Params.dr;
vtkm::FloatDefault Zc = Params.eq_min_z + (vtkm::FloatDefault)(z_i)*Params.dz;
auto Rc_1 = Rc + Params.dr;
auto Zc_1 = Zc + Params.dz;
Rc = (Rc + Rc_1) * 0.5;
Zc = (Zc + Zc_1) * 0.5;
//Get the coeffcients (z,r,4,4)
vtkm::Matrix<vtkm::FloatDefault, 4, 4> acoeff;
//offset = ri * nz + zi
//vtkm::Id offset = (r_i * ncoeff + z_i) * 16; //DRP
vtkm::Id offset = (z_i * Params.ncoeff_r + r_i) * 16;
vtkm::FloatDefault psi, dpsi_dr, dpsi_dz, d2psi_d2r, d2psi_drdz, d2psi_d2z;
EvalBicub2(
R, Z, Rc, Zc, offset, Coeff_2D, psi, dpsi_dr, dpsi_dz, d2psi_drdz, d2psi_d2r, d2psi_d2z);
if (psi < 0)
psi = 0;
metadata.Psi = psi;
//PSI = psi;
metadata.gradPsi_rzp[0] = dpsi_dr;
metadata.gradPsi_rzp[1] = dpsi_dz;
metadata.gradPsi_rzp[2] = 0;
metadata.dpsi_dr = dpsi_dr;
metadata.dpsi_dz = dpsi_dz;
metadata.d2psi_drdz = d2psi_drdz;
metadata.d2psi_d2r = d2psi_d2r;
metadata.d2psi_d2z = d2psi_d2z;
vtkm::FloatDefault fld_I, fld_dIdpsi;
bool isR12 = IsRegion12(R, Z, psi, Params);
if (!isR12)
{
fld_I = I_interpol(psi, 0, 3, Params, Coeff_1D);
fld_dIdpsi = I_interpol(psi, 1, 3, Params, Coeff_1D);
}
else
{
fld_I = I_interpol(psi, 0, 1, Params, Coeff_1D);
fld_dIdpsi = I_interpol(psi, 1, 1, Params, Coeff_1D);
}
vtkm::FloatDefault over_r = 1 / R;
vtkm::FloatDefault over_r2 = over_r * over_r;
vtkm::FloatDefault Br = -dpsi_dz * over_r;
vtkm::FloatDefault Bz = dpsi_dr * over_r;
vtkm::FloatDefault Bp = fld_I * over_r;
metadata.B0_rzp = vtkm::Vec3f(Br, Bz, Bp);
const vtkm::FloatDefault bp_sign = 1.0;
auto dBr_dr = (dpsi_dz * over_r2 - d2psi_drdz * over_r) * bp_sign;
auto dBr_dz = -d2psi_d2z * over_r * bp_sign;
auto dBr_dp = 0.0 * bp_sign;
auto dBz_dr = (-dpsi_dr * over_r2 + d2psi_d2r * over_r) * bp_sign;
auto dBz_dz = d2psi_drdz * over_r * bp_sign;
auto dBz_dp = 0.0 * bp_sign;
auto dBp_dr = dpsi_dr * fld_dIdpsi * over_r - fld_I * over_r2;
auto dBp_dz = fld_dIdpsi * dpsi_dz * over_r;
//auto dBp_dp = jacobian_rzp[2][2];
//calculate curl_B
//vtkm::Vec3f curlB_rzp;
metadata.curlB_rzp[0] = dBz_dp * over_r - dBp_dz;
metadata.curlB_rzp[1] = Bp * over_r + dBp_dr - dBr_dp * over_r;
metadata.curlB_rzp[2] = dBr_dz - dBz_dr;
//std::cout<<"curl_B_rzp= "<<curlB_rzp<<std::endl;
//calculate curl_nb
vtkm::FloatDefault Bmag = vtkm::Magnitude(metadata.B0_rzp);
vtkm::FloatDefault over_B = 1. / Bmag, over_B2 = over_B * over_B;
vtkm::FloatDefault dBdr, dBdz;
dBdr = (Br * dBr_dr + Bp * dBp_dr + Bz * dBz_dr) * over_B;
dBdz = (Br * dBr_dz + Bp * dBp_dz + Bz * dBz_dz) * over_B;
//vtkm::Vec3f curl_nb_rzp;
metadata.curl_nb_rzp[0] = metadata.curlB_rzp[0] * over_B + (Bp * dBdz) * over_B2;
metadata.curl_nb_rzp[1] = metadata.curlB_rzp[1] * over_B + (-Bp * dBdr) * over_B2;
metadata.curl_nb_rzp[2] = metadata.curlB_rzp[2] * over_B + (Bz * dBdr - Br * dBdz) * over_B2;
return true;
}
template <typename ScalarPortal>
VTKM_EXEC vtkm::Vec3f HighOrderBOnly(const vtkm::Vec3f& ptRPZ,
const XGCParams& Params,
const ScalarPortal& Coeff_1D,
const ScalarPortal& Coeff_2D)
{
vtkm::FloatDefault R = ptRPZ[0], Z = ptRPZ[2];
vtkm::Vec3f ptRZ(R, Z, 0);
int r_i = GetIndex(R, Params.nr, Params.rmin, Params.dr_inv);
int z_i = GetIndex(Z, Params.nz, Params.zmin, Params.dz_inv);
// rc(i), zc(j)
vtkm::FloatDefault Rc = Params.rmin + (vtkm::FloatDefault)(r_i)*Params.dr;
vtkm::FloatDefault Zc = Params.zmin + (vtkm::FloatDefault)(z_i)*Params.dz;
auto Rc_1 = Rc + Params.dr;
auto Zc_1 = Zc + Params.dz;
Rc = (Rc + Rc_1) * 0.5;
Zc = (Zc + Zc_1) * 0.5;
//Get the coeffcients (z,r,4,4)
vtkm::Matrix<vtkm::FloatDefault, 4, 4> acoeff;
//vtkm::Id offset = (r_i * ncoeff + z_i) * 16; //DRP
vtkm::Id offset = (z_i * Params.ncoeff_r + r_i) * 16;
/*
std::cout<<"InterpolatePsi: "<<vtkm::Vec2f(R,Z)<<std::endl;
std::cout<<" i/j= "<<r_i<<" "<<z_i<<std::endl;
std::cout<<" ncoeff= "<<ncoeff<<std::endl;
std::cout<<" offset= "<<offset<<std::endl;
std::cout<<" Rc/Zc= "<<Rc<<" "<<Zc<<std::endl;
*/
vtkm::FloatDefault psi, dpsi_dr, dpsi_dz, d2psi_d2r, d2psi_drdz, d2psi_d2z;
EvalBicub2(
R, Z, Rc, Zc, offset, Coeff_2D, psi, dpsi_dr, dpsi_dz, d2psi_drdz, d2psi_d2r, d2psi_d2z);
if (psi < 0)
psi = 0;
bool isR12 = IsRegion12(R, Z, psi, Params);
vtkm::FloatDefault fld_I;
if (!isR12)
fld_I = I_interpol(psi, 0, 3, Params, Coeff_1D);
else
fld_I = I_interpol(psi, 0, 0, Params, Coeff_1D);
vtkm::FloatDefault over_r = 1 / R;
vtkm::FloatDefault Br = -dpsi_dz * over_r;
vtkm::FloatDefault Bz = dpsi_dr * over_r;
vtkm::FloatDefault Bp = fld_I * over_r;
return vtkm::Vec3f(Br, Bz, Bp);
}
template <typename ScalarPortal>
VTKM_EXEC vtkm::Vec3f GetB(vtkm::Vec3f& pt_rpz,
const XGCParams& Params,
const ScalarPortal& Coeff_1D,
const ScalarPortal& Coeff_2D)
{
auto B0_rzp = HighOrderBOnly(pt_rpz, Params, Coeff_1D, Coeff_2D);
vtkm::Vec3f B0_rpz(B0_rzp[0], B0_rzp[2], B0_rzp[1]);
// This gives is the time derivative: Br = dR/dt, Bz= dZ/dt, B_phi/R = dphi/dt
// We need with respect to phi:
// dR/dphi = dR/dt / (dphi/dt) = Br / B_phi * R
// same for z;
B0_rpz[0] /= B0_rzp[2];
B0_rpz[2] /= B0_rzp[2];
B0_rpz[0] *= pt_rpz[0];
B0_rpz[2] *= pt_rpz[0];
return B0_rpz;
}
template <typename ScalarPortal>
VTKM_EXEC bool CalcFieldFollowingPt(const vtkm::Vec3f& pt_rpz,
const vtkm::Vec3f& B0_rpz,
const vtkm::FloatDefault& Phi0,
const vtkm::FloatDefault& Phi1,
const XGCParams& Params,
const ScalarPortal& coeff_1D,
const ScalarPortal& coeff_2D,
vtkm::Vec3f& x_ff_rpz)
{
using Ray3f = vtkm::Ray<vtkm::FloatDefault, 3, true>;
vtkm::FloatDefault R = pt_rpz[0];
vtkm::FloatDefault Phi = pt_rpz[1];
vtkm::FloatDefault Z = pt_rpz[2];
vtkm::FloatDefault PhiMid = Phi0 + (Phi1 - Phi0) / 2.0;
vtkm::Plane<> midPlane({ 0, PhiMid, 0 }, { 0, 1, 0 });
Ray3f ray_rpz({ R, Phi, Z }, B0_rpz);
//Get point on mid plane. Use the R,Z for this point for triangle finds.
vtkm::FloatDefault RP_T;
vtkm::Vec3f ptOnMidPlane_rpz;
bool b;
midPlane.Intersect(ray_rpz, RP_T, ptOnMidPlane_rpz, b);
//Now, do it using RK4 and two steps.
vtkm::FloatDefault h = (PhiMid - Phi) / 2.0;
//h = -h;
vtkm::FloatDefault h_2 = h / 2.0;
//k1 = F(p)
//k2 = F(p+hk1/2)
//k3 = F(p+hk2/2)
//k4 = F(p+hk3)
//Yn+1 = Yn + 1/6 h (k1+2k2+2k3+k4)
vtkm::Vec3f p0 = { R, Phi, Z }; //pt_rpz;
vtkm::Vec3f tmp, k1, k2, k3, k4;
//std::cout<<" p0 = "<<p0<<std::endl;
for (int i = 0; i < 2; i++)
{
k1 = GetB(p0, Params, coeff_1D, coeff_2D);
tmp = p0 + k1 * h_2;
k2 = GetB(tmp, Params, coeff_1D, coeff_2D);
tmp = p0 + k2 * h_2;
k3 = GetB(tmp, Params, coeff_1D, coeff_2D);
tmp = p0 + k3 * h;
k4 = GetB(tmp, Params, coeff_1D, coeff_2D);
vtkm::Vec3f vec = (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;
p0 = p0 + h * vec;
}
x_ff_rpz = p0;
return true;
}
template <typename PortalType>
VTKM_EXEC typename PortalType::ValueType Eval(const PortalType& portal,
const vtkm::Id& offset,
const vtkm::Vec3f& param,
const vtkm::Vec<vtkm::Id, 3>& vId)
{
#if 1
//Hard code...
const auto& v0 = portal.Get(vId[0] + offset);
const auto& v1 = portal.Get(vId[1] + offset);
const auto& v2 = portal.Get(vId[2] + offset);
const auto& w1 = param[0];
const auto& w2 = param[1];
const auto w0 = 1 - (w1 + w2);
auto ret = v0 * w0 + v1 * w1 + v2 * w2;
#else
T ret;
vtkm::VecVariable<T, 3> vals;
vals.Append(portal.Get(vId[0] + offset));
vals.Append(portal.Get(vId[1] + offset));
vals.Append(portal.Get(vId[2] + offset));
vtkm::exec::CellInterpolate(vals, param, vtkm::CellShapeTagTriangle(), ret);
#endif
return ret;
}
VTKM_EXEC inline void GetPlaneIdx(const vtkm::FloatDefault& phi,
const vtkm::FloatDefault& Period,
const XGCParams& Params,
vtkm::FloatDefault& phiN,
vtkm::Id& plane0,
vtkm::Id& plane1,
vtkm::FloatDefault& phi0,
vtkm::FloatDefault& phi1,
vtkm::Id& numRevs,
vtkm::FloatDefault& T)
{
numRevs = vtkm::Floor(vtkm::Abs(phi / Period));
phiN = phi;
if (phi < 0)
{
phiN += ((1 + numRevs) * Period);
}
else if (phi > Period)
{
phiN -= (numRevs * Period);
}
plane0 = vtkm::Floor(phiN / Params.dPhi);
if (plane0 == Params.NumPlanes)
plane0 = 0;
plane1 = plane0 + 1;
phi0 = static_cast<vtkm::FloatDefault>(plane0) * Params.dPhi;
phi1 = static_cast<vtkm::FloatDefault>(plane1) * Params.dPhi;
if (plane1 == Params.NumPlanes)
plane1 = 0;
T = (phiN - phi0) / (phi1 - phi0);
}
}
}
} //vtkm::worklet::flow
#endif //vtkm_filter_flow_worklet_XGCHelper_h

284
vtkm/filter/flow/worklet/XGCPoincare.h Normal file
View File

@ -0,0 +1,284 @@
//=============================================================================
//
// 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 vtkm_worklet_particleadvection_xgcpoincare_h
#define vtkm_worklet_particleadvection_xgcpoincare_h
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/filter/flow/worklet/XGCHelper.h>
namespace vtkm
{
namespace worklet
{
namespace flow
{
template <typename ParticleType>
class XGCPoincareExec
{
public:
VTKM_EXEC_CONT
XGCPoincareExec()
: NumParticles(0)
, MaxSteps(0)
, MaxPunctures(0)
, PunctureCounts()
, Validity()
, PunctureIDs()
, OutputR()
, OutputZ()
, OutputTheta()
, OutputPsi()
{
}
VTKM_CONT
XGCPoincareExec(vtkm::Id numParticles,
vtkm::Id maxSteps,
vtkm::Id maxPunctures,
const XGCParams& params,
vtkm::FloatDefault period,
const vtkm::cont::ArrayHandle<vtkm::FloatDefault>& coeff_1D,
const vtkm::cont::ArrayHandle<vtkm::FloatDefault>& coeff_2D,
const vtkm::cont::ArrayHandle<vtkm::Id>& punctureCounts,
const vtkm::cont::ArrayHandle<vtkm::Id>& validity,
const vtkm::cont::ArrayHandle<vtkm::Id>& punctureIDs,
const vtkm::cont::ArrayHandle<vtkm::FloatDefault>& outputR,
const vtkm::cont::ArrayHandle<vtkm::FloatDefault>& outputZ,
const vtkm::cont::ArrayHandle<vtkm::FloatDefault>& outputTheta,
const vtkm::cont::ArrayHandle<vtkm::FloatDefault>& outputPsi,
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token)
: NumParticles(numParticles)
, MaxSteps(maxSteps)
, MaxPunctures(maxPunctures)
, Params(params)
, Period(period)
{
Coeff_1D = coeff_1D.PrepareForInput(device, token);
Coeff_2D = coeff_2D.PrepareForInput(device, token);
PunctureCounts = punctureCounts.PrepareForInPlace(device, token);
Validity = validity.PrepareForInPlace(device, token);
PunctureIDs = punctureIDs.PrepareForInPlace(device, token);
OutputR = outputR.PrepareForOutput(this->NumParticles * this->MaxPunctures, device, token);
OutputZ = outputZ.PrepareForOutput(this->NumParticles * this->MaxPunctures, device, token);
OutputTheta =
outputTheta.PrepareForOutput(this->NumParticles * this->MaxPunctures, device, token);
OutputPsi = outputPsi.PrepareForOutput(this->NumParticles * this->MaxPunctures, device, token);
}
//template <typename ParticleType>
VTKM_EXEC void Analyze(const vtkm::Id index,
const ParticleType& oldParticle,
ParticleType& newParticle) const
{
vtkm::Id numRevs0 = vtkm::Floor(vtkm::Abs(oldParticle.Pos[1] / this->Period));
vtkm::Id numRevs1 = vtkm::Floor(vtkm::Abs(newParticle.Pos[1] / this->Period));
//if (this->SaveTraces)
// Intersections.Set(idx*this->MaxIter + particle.NumSteps, particle.Pos);
//std::cout<<std::setprecision(12)<<" Step: "<<particle.Pos<<std::endl;
if (numRevs1 > numRevs0)
{
std::cout << "Puncture!!! : " << newParticle.Pos << std::endl;
auto R = newParticle.Pos[0];
auto Z = newParticle.Pos[2];
auto theta = vtkm::ATan2(Z - this->Params.eq_axis_z, R - this->Params.eq_axis_r);
if (theta < 0)
theta += vtkm::TwoPi();
//calcualte psi. need to stash psi on the particle somehow....
vtkm::Vec3f ptRPZ = newParticle.Pos;
ParticleMetaData metadata;
HighOrderB(ptRPZ, metadata, this->Params, this->Coeff_1D, this->Coeff_2D);
vtkm::Id loc = (index * this->MaxPunctures) + newParticle.NumPunctures;
std::cout << "Puncture data : " << R << " " << Z << " " << theta << " "
<< (metadata.Psi / this->Params.eq_x_psi) << " " << index << std::endl;
this->Validity.Set(loc, 1);
this->OutputR.Set(loc, R);
this->OutputZ.Set(loc, Z);
this->OutputTheta.Set(loc, theta);
this->OutputPsi.Set(loc, metadata.Psi / this->Params.eq_x_psi);
this->PunctureIDs.Set(loc, index);
newParticle.NumPunctures++;
}
}
private:
using ScalarPortalR = typename vtkm::cont::ArrayHandle<vtkm::FloatDefault>::ReadPortalType;
using ScalarPortalW = typename vtkm::cont::ArrayHandle<vtkm::FloatDefault>::WritePortalType;
using IdPortal = typename vtkm::cont::ArrayHandle<vtkm::Id>::WritePortalType;
using VecPortal = typename vtkm::cont::ArrayHandle<vtkm::Vec3f>::WritePortalType;
vtkm::Id NumParticles;
vtkm::Id MaxSteps;
vtkm::Id MaxPunctures;
vtkm::FloatDefault Period;
XGCParams Params;
ScalarPortalR Coeff_1D;
ScalarPortalR Coeff_2D;
ScalarPortalW OutputR;
ScalarPortalW OutputZ;
ScalarPortalW OutputTheta;
ScalarPortalW OutputPsi;
IdPortal PunctureIDs;
//VecPortal Intersections;
IdPortal PunctureCounts;
IdPortal Validity;
};
template <typename ParticleType>
class XGCPoincare : public vtkm::cont::ExecutionObjectBase
{
public:
// Intended to store advected particles after Finalize
vtkm::cont::ArrayHandle<ParticleType> Particles;
vtkm::cont::ArrayHandle<vtkm::Id> PunctureIDs;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> OutputR;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> OutputZ;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> OutputTheta;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> OutputPsi;
//Helper functor for compacting history
struct IsOne
{
template <typename T>
VTKM_EXEC_CONT bool operator()(const T& x) const
{
return x == T(1);
}
};
VTKM_CONT
XGCPoincare(vtkm::Id maxSteps,
vtkm::Id maxPunctures,
const vtkm::cont::ArrayHandle<vtkm::FloatDefault>& coeff_1D,
const vtkm::cont::ArrayHandle<vtkm::FloatDefault>& coeff_2D,
const XGCParams& params,
vtkm::FloatDefault period)
: Particles()
, MaxSteps(maxSteps)
, MaxPunctures(maxPunctures)
, Coeff_1D(coeff_1D)
, Coeff_2D(coeff_2D)
, Params(params)
, Period(period)
{
}
VTKM_CONT
void UseAsTemplate(const XGCPoincare& other)
{
this->MaxSteps = other.MaxSteps;
this->MaxPunctures = other.MaxPunctures;
this->Period = other.Period;
this->Params = other.Params;
this->Coeff_1D = other.Coeff_1D;
this->Coeff_2D = other.Coeff_2D;
}
VTKM_CONT XGCPoincareExec<ParticleType> PrepareForExecution(vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
{
return XGCPoincareExec<ParticleType>(this->NumParticles,
this->MaxSteps,
this->MaxPunctures,
this->Params,
this->Period,
this->Coeff_1D,
this->Coeff_2D,
this->PunctureCounts,
this->Validity,
this->PunctureIDs,
this->OutputR,
this->OutputZ,
this->OutputTheta,
this->OutputPsi,
device,
token);
}
VTKM_CONT bool SupportPushOutOfBounds() const { return false; }
VTKM_CONT
void InitializeAnalysis(const vtkm::cont::ArrayHandle<ParticleType>& particles)
{
this->NumParticles = particles.GetNumberOfValues();
//Create Validity initialized to zero.
vtkm::cont::ArrayHandleConstant<vtkm::Id> validity(0, this->NumParticles * this->MaxPunctures);
vtkm::cont::ArrayCopy(validity, this->Validity);
//Create PunctureCount initialized to zero.
vtkm::cont::ArrayHandleConstant<vtkm::Id> punctureCounts(0, this->NumParticles);
vtkm::cont::ArrayCopy(punctureCounts, this->PunctureCounts);
vtkm::cont::ArrayHandleConstant<vtkm::Id> initIds(-1, this->NumParticles * this->MaxPunctures);
vtkm::cont::ArrayCopy(initIds, this->PunctureIDs);
}
VTKM_CONT
void FinalizeAnalysis(vtkm::cont::ArrayHandle<ParticleType>& particles)
{
(void)particles;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> outR, outZ, outTheta, outPsi;
vtkm::cont::ArrayHandle<vtkm::Id> outIds;
vtkm::cont::Algorithm::CopyIf(this->OutputR, this->Validity, outR, IsOne());
vtkm::cont::Algorithm::CopyIf(this->OutputZ, this->Validity, outZ, IsOne());
vtkm::cont::Algorithm::CopyIf(this->OutputTheta, this->Validity, outTheta, IsOne());
vtkm::cont::Algorithm::CopyIf(this->OutputPsi, this->Validity, outPsi, IsOne());
vtkm::cont::Algorithm::CopyIf(this->PunctureIDs, this->Validity, outIds, IsOne());
vtkm::cont::ArrayCopyShallowIfPossible(outR, this->OutputR);
vtkm::cont::ArrayCopyShallowIfPossible(outZ, this->OutputZ);
vtkm::cont::ArrayCopyShallowIfPossible(outTheta, this->OutputTheta);
vtkm::cont::ArrayCopyShallowIfPossible(outPsi, this->OutputPsi);
vtkm::cont::ArrayCopyShallowIfPossible(outIds, this->PunctureIDs);
}
VTKM_CONT static bool MakeDataSet(vtkm::cont::DataSet& dataset,
const std::vector<XGCPoincare>& results)
{
(void)dataset;
(void)results;
return true;
}
private:
vtkm::Id NumParticles;
vtkm::Id MaxSteps;
vtkm::Id MaxPunctures;
// For Higher order interpolation
// on RZ uniform grid, not the triangle grid.
vtkm::cont::ArrayHandle<vtkm::FloatDefault> Coeff_1D;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> Coeff_2D;
XGCParams Params;
vtkm::FloatDefault Period;
// Output Arrays
vtkm::cont::ArrayHandle<vtkm::Id> PunctureCounts;
//vtkm::cont::ArrayHandle<vtkm::Id> InitialPunctureCounts;
vtkm::cont::ArrayHandle<vtkm::Id> Validity;
};
} // namespace particleadvection
} // namespace worklet
} // namespace vtkm
#endif //vtkm_worklet_particleadvection_xgcpoincare_h