Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-16 17:22:55 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Unit tests for AOS particle advection.

Browse Source
This commit is contained in:
Dave Pugmire 2019-10-01 13:33:12 -04:00
parent 28a73494b1
commit 6e9ab1ce5d
9 changed files with 343 additions and 237 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
vtkm/Bitset.h
View File

@ -41,7 +41,7 @@ struct Bitset
this->Mask = this->Mask ^ (static_cast<MaskType>(0) << bitIndex);
}
VTKM_EXEC_CONT bool test(vtkm::Id bitIndex)
VTKM_EXEC_CONT bool test(vtkm::Id bitIndex) const
{
return ((this->Mask & (static_cast<MaskType>(1) << bitIndex)) != 0);
}

66
vtkm/Particle.h
View File

@ -10,22 +10,58 @@
#ifndef vtk_m_Particle_h
#define vtk_m_Particle_h
#include <vtkm/Bitset.h>
namespace vtkm
{
enum ParticleStatus : vtkm::Id
//Bit field describing the status:
class ParticleStatus : public vtkm::Bitset<vtkm::UInt8>
{
UNDEFINED = 0,
SUCCESS = 1,
TERMINATED = 1 << 1,
EXIT_SPATIAL_BOUNDARY = 1 << 2,
EXIT_TEMPORAL_BOUNDARY = 1 << 3,
FAIL = 1 << 4,
TOOK_ANY_STEPS = 1 << 5
public:
VTKM_EXEC_CONT ParticleStatus()
{
this->SetOk();
this->ClearTerminate();
}
VTKM_EXEC_CONT void SetOk() { this->set(SUCCESS_BIT); }
VTKM_EXEC_CONT bool CheckOk() const { return this->test(SUCCESS_BIT); }
VTKM_EXEC_CONT void SetFail() { this->reset(SUCCESS_BIT); }
VTKM_EXEC_CONT bool CheckFail() const { return !this->test(SUCCESS_BIT); }
VTKM_EXEC_CONT void SetTerminate() { this->set(TERMINATE_BIT); }
VTKM_EXEC_CONT void ClearTerminate() { this->reset(TERMINATE_BIT); }
VTKM_EXEC_CONT bool CheckTerminate() const { return this->test(TERMINATE_BIT); }
VTKM_EXEC_CONT void SetSpatialBounds() { this->set(SPATIAL_BOUNDS_BIT); }
VTKM_EXEC_CONT void ClearSpatialBounds() { this->reset(SPATIAL_BOUNDS_BIT); }
VTKM_EXEC_CONT bool CheckSpatialBounds() const { return this->test(SPATIAL_BOUNDS_BIT); }
VTKM_EXEC_CONT void SetTemporalBounds() { this->set(TEMPORAL_BOUNDS_BIT); }
VTKM_EXEC_CONT void ClearTemporalBounds() { this->reset(TEMPORAL_BOUNDS_BIT); }
VTKM_EXEC_CONT bool CheckTemporalBounds() const { return this->test(TEMPORAL_BOUNDS_BIT); }
VTKM_EXEC_CONT void SetTookAnySteps() { this->set(TOOK_ANY_STEPS_BIT); }
VTKM_EXEC_CONT void ClearTookAnySteps() { this->reset(TOOK_ANY_STEPS_BIT); }
VTKM_EXEC_CONT bool CheckTookAnySteps() const { return this->test(TOOK_ANY_STEPS_BIT); }
private:
static constexpr vtkm::IdComponent SUCCESS_BIT = 0;
static constexpr vtkm::IdComponent TERMINATE_BIT = 1;
static constexpr vtkm::IdComponent SPATIAL_BOUNDS_BIT = 2;
static constexpr vtkm::IdComponent TEMPORAL_BOUNDS_BIT = 3;
static constexpr vtkm::IdComponent TOOK_ANY_STEPS_BIT = 4;
};
inline VTKM_EXEC_CONT std::ostream& operator<<(std::ostream& s, const vtkm::ParticleStatus& status)
{
s << "[" << status.CheckOk() << " " << status.CheckTerminate() << " "
<< status.CheckSpatialBounds() << " " << status.CheckTemporalBounds() << "]";
return s;
}
class Particle
{
public:
@ -34,7 +70,7 @@ public:
: Pos()
, ID(-1)
, NumSteps(0)
, Status(vtkm::ParticleStatus::UNDEFINED)
, Status()
, Time(0)
{
}
@ -51,10 +87,10 @@ public:
VTKM_EXEC_CONT
Particle(const vtkm::Vec3f& p,
vtkm::Id id,
vtkm::Id numSteps,
vtkm::Id status = vtkm::ParticleStatus::UNDEFINED,
vtkm::FloatDefault time = 0)
const vtkm::Id& id,
const vtkm::Id& numSteps = 0,
const vtkm::ParticleStatus& status = vtkm::ParticleStatus(),
const vtkm::FloatDefault& time = 0)
: Pos(p)
, ID(id)
, NumSteps(numSteps)
@ -66,7 +102,7 @@ public:
vtkm::Vec3f Pos;
vtkm::Id ID;
vtkm::Id NumSteps;
vtkm::Id Status;
vtkm::ParticleStatus Status;
vtkm::FloatDefault Time;
};
}

8
vtkm/worklet/ParticleAdvection.h
View File

@ -139,15 +139,11 @@ public:
template <typename IntegratorType, typename ParticleStorage>
ParticleAdvectionResult Run(const IntegratorType& it,
vtkm::cont::ArrayHandle<vtkm::Particle, ParticleStorage>& particles,
const vtkm::Id& MaxSteps,
int which)
const vtkm::Id& MaxSteps)
{
vtkm::worklet::particleadvection::ParticleAdvectionWorklet<IntegratorType> worklet;
if (which == 0)
worklet.Run(it, particles, MaxSteps);
else
worklet.Run2(it, particles, MaxSteps);
worklet.Run(it, particles, MaxSteps);
return ParticleAdvectionResult();
}

1
vtkm/worklet/particleadvection/CMakeLists.txt
View File

@ -10,7 +10,6 @@
set(headers
CellInterpolationHelper.h
EvaluatorStatus.h
GridEvaluators.h
Integrators.h
Particles.h

8
vtkm/worklet/particleadvection/EvaluatorStatus.h
View File

@ -13,18 +13,14 @@
#ifndef vtk_m_worklet_particleadvection_EvaluatorStatus_h
#define vtk_m_worklet_particleadvection_EvaluatorStatus_h
#include <vtkm/Bitset.h>
namespace vtkm
{
namespace worklet
{
namespace particleadvection
{
enum class EvaluatorStatus
{
SUCCESS = 0,
OUTSIDE_SPATIAL_BOUNDS,
OUTSIDE_TEMPORAL_BOUNDS
};
} //namespace particleadvection
} //namespace worklet

33
vtkm/worklet/particleadvection/GridEvaluators.h
View File

@ -11,6 +11,7 @@
#ifndef vtk_m_worklet_particleadvection_GridEvaluators_h
#define vtk_m_worklet_particleadvection_GridEvaluators_h
#include <vtkm/Bitset.h>
#include <vtkm/Types.h>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/cont/ArrayHandle.h>
@ -23,7 +24,6 @@
#include <vtkm/cont/DeviceAdapter.h>
#include <vtkm/worklet/particleadvection/CellInterpolationHelper.h>
#include <vtkm/worklet/particleadvection/EvaluatorStatus.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
namespace vtkm
@ -32,6 +32,26 @@ namespace worklet
{
namespace particleadvection
{
class EvaluatorStatus : public vtkm::Bitset<vtkm::UInt8>
{
public:
void SetOk() { this->set(SUCCESS_BIT); }
bool CheckOk() const { return this->test(SUCCESS_BIT); }
void SetFail() { this->reset(SUCCESS_BIT); }
bool CheckFail() const { return !this->test(SUCCESS_BIT); }
void SetSpatialBounds() { this->set(SPATIAL_BOUNDS_BIT); }
bool CheckSpatialBounds() const { return this->test(SPATIAL_BOUNDS_BIT); }
void SetTemporalBounds() { this->set(TEMPORAL_BOUNDS_BIT); }
bool CheckTemporalBounds() const { return this->test(TEMPORAL_BOUNDS_BIT); }
private:
static constexpr vtkm::IdComponent SUCCESS_BIT = 0;
static constexpr vtkm::IdComponent SPATIAL_BOUNDS_BIT = 1;
static constexpr vtkm::IdComponent TEMPORAL_BOUNDS_BIT = 2;
};
template <typename DeviceAdapter, typename FieldArrayType>
class ExecutionGridEvaluator
@ -94,9 +114,15 @@ public:
vtkm::Id cellId;
Point parametric;
vtkm::exec::FunctorBase tmp;
EvaluatorStatus status;
Locator->FindCell(point, cellId, parametric, tmp);
if (cellId == -1)
return EvaluatorStatus::OUTSIDE_SPATIAL_BOUNDS;
{
status.SetFail();
status.SetSpatialBounds();
return status;
}
vtkm::UInt8 cellShape;
vtkm::IdComponent nVerts;
@ -108,7 +134,8 @@ public:
fieldValues.Append(Field.Get(ptIndices[i]));
out = vtkm::exec::CellInterpolate(fieldValues, parametric, cellShape, tmp);
return EvaluatorStatus::SUCCESS;
status.SetOk();
return status;
}
private:

206
vtkm/worklet/particleadvection/Integrators.h
View File

@ -13,8 +13,10 @@
#ifndef vtk_m_worklet_particleadvection_Integrators_h
#define vtk_m_worklet_particleadvection_Integrators_h
#include <iomanip>
#include <limits>
#include <vtkm/Bitset.h>
#include <vtkm/TypeTraits.h>
#include <vtkm/Types.h>
#include <vtkm/VectorAnalysis.h>
@ -22,7 +24,7 @@
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/VirtualObjectHandle.h>
#include <vtkm/worklet/particleadvection/EvaluatorStatus.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Particles.h>
namespace vtkm
@ -31,14 +33,58 @@ namespace worklet
{
namespace particleadvection
{
enum class IntegratorStatus
class IntegratorStatus : public vtkm::Bitset<vtkm::UInt8>
{
SUCCESS = 0,
OUTSIDE_SPATIAL_BOUNDS,
OUTSIDE_TEMPORAL_BOUNDS,
FAIL
public:
VTKM_EXEC_CONT IntegratorStatus() {}
VTKM_EXEC_CONT IntegratorStatus(const bool& ok, const bool& spatial, const bool& temporal)
{
SetBit(0 /*SUCCESS_BIT*/, ok);
SetBit(1 /*SPATIAL_BOUNDS_BIT*/, spatial);
SetBit(2 /*TEMPORAL_BOUNDS_BIT*/, temporal);
}
VTKM_EXEC_CONT IntegratorStatus(const EvaluatorStatus& es)
{
SetBit(0 /*SUCCESS_BIT*/, es.CheckOk());
SetBit(1 /*SPATIAL_BOUNDS_BIT*/, es.CheckSpatialBounds());
SetBit(2 /*TEMPORAL_BOUNDS_BIT*/, es.CheckTemporalBounds());
}
VTKM_EXEC_CONT void SetOk() { this->set(SUCCESS_BIT); }
VTKM_EXEC_CONT bool CheckOk() const { return this->test(SUCCESS_BIT); }
VTKM_EXEC_CONT void SetFail() { this->reset(SUCCESS_BIT); }
VTKM_EXEC_CONT bool CheckFail() const { return !this->test(SUCCESS_BIT); }
VTKM_EXEC_CONT void SetSpatialBounds() { this->set(SPATIAL_BOUNDS_BIT); }
VTKM_EXEC_CONT bool CheckSpatialBounds() const { return this->test(SPATIAL_BOUNDS_BIT); }
VTKM_EXEC_CONT void SetTemporalBounds() { this->set(TEMPORAL_BOUNDS_BIT); }
VTKM_EXEC_CONT bool CheckTemporalBounds() const { return this->test(TEMPORAL_BOUNDS_BIT); }
private:
VTKM_EXEC_CONT void SetBit(const vtkm::IdComponent& bit, bool val)
{
if (val)
this->set(bit);
else
this->reset(bit);
}
const static vtkm::IdComponent SUCCESS_BIT = 0;
const static vtkm::IdComponent SPATIAL_BOUNDS_BIT = 1;
const static vtkm::IdComponent TEMPORAL_BOUNDS_BIT = 2;
};
inline VTKM_EXEC_CONT std::ostream& operator<<(std::ostream& s, const IntegratorStatus& status)
{
s << "[" << status.CheckOk() << " " << status.CheckSpatialBounds() << " "
<< status.CheckTemporalBounds() << "]";
return s;
}
class Integrator : public vtkm::cont::ExecutionObjectBase
{
protected:
@ -73,22 +119,6 @@ public:
vtkm::FloatDefault& time,
vtkm::Vec3f& outpos) const = 0;
VTKM_EXEC
IntegratorStatus ConvertToIntegratorStatus(EvaluatorStatus status) const
{
switch (status)
{
case EvaluatorStatus::SUCCESS:
return IntegratorStatus::SUCCESS;
case EvaluatorStatus::OUTSIDE_SPATIAL_BOUNDS:
return IntegratorStatus::OUTSIDE_SPATIAL_BOUNDS;
case EvaluatorStatus::OUTSIDE_TEMPORAL_BOUNDS:
return IntegratorStatus::OUTSIDE_TEMPORAL_BOUNDS;
default:
return IntegratorStatus::FAIL;
}
}
protected:
vtkm::FloatDefault StepLength = 1.0f;
vtkm::FloatDefault Tolerance = 0.001f;
@ -135,14 +165,21 @@ protected:
{
// If particle is out of either spatial or temporal boundary to begin with,
// then return the corresponding status.
IntegratorStatus status;
if (!this->Evaluator.IsWithinSpatialBoundary(inpos))
return IntegratorStatus::OUTSIDE_SPATIAL_BOUNDS;
{
status.SetSpatialBounds();
return status;
}
if (!this->Evaluator.IsWithinTemporalBoundary(time))
return IntegratorStatus::OUTSIDE_TEMPORAL_BOUNDS;
{
status.SetTemporalBounds();
return status;
}
vtkm::Vec3f velocity;
IntegratorStatus status = CheckStep(inpos, this->StepLength, time, velocity);
if (status == IntegratorStatus::SUCCESS)
status = CheckStep(inpos, this->StepLength, time, velocity);
if (status.CheckOk())
{
outpos = inpos + StepLength * velocity;
time += StepLength;
@ -159,54 +196,55 @@ protected:
vtkm::Vec3f& outpos) const override
{
if (!this->Evaluator.IsWithinSpatialBoundary(inpos))
return IntegratorStatus::OUTSIDE_SPATIAL_BOUNDS;
return IntegratorStatus(false, true, false);
if (!this->Evaluator.IsWithinTemporalBoundary(time))
return IntegratorStatus::OUTSIDE_TEMPORAL_BOUNDS;
return IntegratorStatus(false, false, true);
vtkm::FloatDefault optimalLength = static_cast<vtkm::FloatDefault>(0);
vtkm::Id iteration = static_cast<vtkm::Id>(1);
vtkm::Id maxIterations = static_cast<vtkm::Id>(1 << 20);
vtkm::Vec3f velocity;
vtkm::Vec3f workpos(inpos);
vtkm::FloatDefault worktime = time;
// According to the earlier checks this call to Evaluate should return
// the correct velocity at the current location, this is to use just in
// case we are not able to find the optimal lenght in 20 iterations..
this->Evaluator.Evaluate(workpos, time, velocity);
while (iteration < maxIterations)
{
iteration = iteration << 1;
vtkm::FloatDefault length =
optimalLength + (this->StepLength / static_cast<vtkm::FloatDefault>(iteration));
IntegratorStatus status = this->CheckStep(inpos, length, time, velocity);
if (status == IntegratorStatus::SUCCESS &&
this->Evaluator.IsWithinSpatialBoundary(inpos + velocity * length))
{
workpos = inpos + velocity * length;
worktime = time + length;
optimalLength = length;
}
}
this->Evaluator.Evaluate(workpos, worktime, velocity);
// We have calculated a large enough step length to push the particle
// using the higher order evaluator, take a step using that evaluator.
// Take one final step, which should be an Euler step just to push the
// particle out of the domain boundary
vtkm::Bounds bounds = this->Evaluator.GetSpatialBoundary();
vtkm::Vec3f direction = velocity / vtkm::Magnitude(velocity);
//Stepping by this->StepLength goes beyond the bounds of the dataset.
//We need to take an Euler step that goes outside of the dataset.
//Use a binary search to find the largest step INSIDE the dataset.
//Binary search uses a shrinking bracket of inside / outside, so when
//we terminate, the outside value is the stepsize that will nudge
//the particle outside the dataset.
//The binary search will be between [0, this->StepLength]
vtkm::FloatDefault stepShort = 0, stepLong = this->StepLength;
vtkm::Vec3f currPos(inpos), velocity;
vtkm::FloatDefault currTime = time;
auto evalStatus = this->Evaluator.Evaluate(currPos, time, velocity);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
const vtkm::FloatDefault eps = vtkm::Epsilon<vtkm::FloatDefault>();
vtkm::FloatDefault xStepLength =
vtkm::Abs(direction[0] * eps * static_cast<vtkm::FloatDefault>(bounds.X.Length()));
vtkm::FloatDefault yStepLength =
vtkm::Abs(direction[1] * eps * static_cast<vtkm::FloatDefault>(bounds.Y.Length()));
vtkm::FloatDefault zStepLength =
vtkm::Abs(direction[2] * eps * static_cast<vtkm::FloatDefault>(bounds.Z.Length()));
vtkm::FloatDefault minLength = vtkm::Min(xStepLength, vtkm::Min(yStepLength, zStepLength));
vtkm::FloatDefault div = 1;
for (int i = 0; i < 50; i++)
{
div *= 2;
vtkm::FloatDefault stepCurr = stepShort + (this->StepLength / div);
//See if we can step by stepCurr.
IntegratorStatus status = this->CheckStep(inpos, stepCurr, time, velocity);
if (status.CheckOk())
{
currPos = inpos + velocity * stepShort;
stepShort = stepCurr;
}
else
stepLong = stepCurr;
outpos = workpos + minLength * velocity;
time = worktime + minLength;
return IntegratorStatus::OUTSIDE_SPATIAL_BOUNDS;
//Stop if step bracket is small enough.
if (stepLong - stepShort < eps)
break;
}
//Take Euler step.
currTime = time + stepShort;
evalStatus = this->Evaluator.Evaluate(currPos, currTime, velocity);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
outpos = currPos + stepLong * velocity;
return IntegratorStatus(true, true, !this->Evaluator.IsWithinTemporalBoundary(time));
}
VTKM_EXEC
@ -299,22 +337,22 @@ public:
vtkm::Vec3f k1 = vtkm::TypeTraits<vtkm::Vec3f>::ZeroInitialization();
vtkm::Vec3f k2 = k1, k3 = k1, k4 = k1;
EvaluatorStatus status;
status = this->Evaluator.Evaluate(inpos, time, k1);
if (status != EvaluatorStatus::SUCCESS)
return this->ConvertToIntegratorStatus(status);
status = this->Evaluator.Evaluate(inpos + var1 * k1, var2, k2);
if (status != EvaluatorStatus::SUCCESS)
return this->ConvertToIntegratorStatus(status);
status = this->Evaluator.Evaluate(inpos + var1 * k2, var2, k3);
if (status != EvaluatorStatus::SUCCESS)
return this->ConvertToIntegratorStatus(status);
status = this->Evaluator.Evaluate(inpos + stepLength * k3, var3, k4);
if (status != EvaluatorStatus::SUCCESS)
return this->ConvertToIntegratorStatus(status);
EvaluatorStatus evalStatus;
evalStatus = this->Evaluator.Evaluate(inpos, time, k1);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
evalStatus = this->Evaluator.Evaluate(inpos + var1 * k1, var2, k2);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
evalStatus = this->Evaluator.Evaluate(inpos + var1 * k2, var2, k3);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
evalStatus = this->Evaluator.Evaluate(inpos + stepLength * k3, var3, k4);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
velocity = (k1 + 2 * k2 + 2 * k3 + k4) / 6.0f;
return IntegratorStatus::SUCCESS;
return IntegratorStatus(true, false, evalStatus.CheckTemporalBounds());
}
};
@ -377,7 +415,7 @@ public:
vtkm::Vec3f& velocity) const
{
EvaluatorStatus status = this->Evaluator.Evaluate(inpos, time, velocity);
return this->ConvertToIntegratorStatus(status);
return IntegratorStatus(status);
}
};

211
vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h
View File

@ -24,7 +24,7 @@
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/Particles.h>
#include <vtkm/worklet/particleadvection/ParticlesAOS.h>
//#include <vtkm/worklet/particleadvection/ParticlesAOS.h>
namespace vtkm
{
@ -42,10 +42,11 @@ public:
using InputDomain = _1;
template <typename IntegratorType, typename ParticleType>
VTKM_EXEC void operator()(const vtkm::Id& idx,
const IntegratorType* integrator,
ParticleType& particle) const
VTKM_EXEC void operator()(const vtkm::Id& /*idx*/,
const IntegratorType* /*integrator*/,
ParticleType& /*particle*/) const
{
#if 0
vtkm::Vec3f inpos = particle.GetPos(idx);
vtkm::Vec3f outpos;
vtkm::FloatDefault time = particle.GetTime(idx);
@ -57,7 +58,7 @@ public:
status = integrator->Step(inpos, time, outpos);
// If the status is OK, we only need to check if the particle
// has completed the maximum steps required.
if (status == IntegratorStatus::SUCCESS)
if (status.CheckOk())
{
particle.TakeStep(idx, outpos, time);
// This is to keep track of the particle's time.
@ -70,26 +71,26 @@ public:
// push it a little out of the boundary so that it will start advection in
// another domain, or in another time slice. Taking small steps enables
// reducing the error introduced at spatial or temporal boundaries.
else if (status == IntegratorStatus::OUTSIDE_TEMPORAL_BOUNDS)
else if (status.CheckTemporalBounds())
{
particle.SetExitTemporalBoundary(idx);
break;
}
else if (status == IntegratorStatus::OUTSIDE_SPATIAL_BOUNDS)
else if (status.CheckSpatialBounds())
{
status = integrator->SmallStep(inpos, time, outpos);
particle.TakeStep(idx, outpos, time);
if (status == IntegratorStatus::OUTSIDE_TEMPORAL_BOUNDS)
if (status.CheckTemporalBounds())
{
particle.SetExitTemporalBoundary(idx);
break;
}
else if (status == IntegratorStatus::OUTSIDE_SPATIAL_BOUNDS)
else if (status.CheckSpatialBounds())
{
particle.SetExitSpatialBoundary(idx);
break;
}
else if (status == IntegratorStatus::FAIL)
else if (status.CheckFail())
{
particle.SetError(idx);
break;
@ -97,19 +98,76 @@ public:
}
}
particle.SetTookAnySteps(idx, tookAnySteps);
#endif
}
};
//for streamlines: inherit from the worklet below.
// overload take step to do something else.
//
class ParticleAdvectWorkletAOS2 : public vtkm::worklet::WorkletMapField
{
public:
using ControlSignature = void(ExecObject integrator, FieldInOut particles, FieldIn maxSteps);
using ExecutionSignature = void(_1, _2, _3);
using InputDomain = _3;
using ControlSignature = void(FieldIn idx,
ExecObject integrator,
FieldInOut particles,
FieldIn maxSteps);
using ExecutionSignature = void(_1 idx, _2 integrator, _3 particles, _4 maxSteps);
using InputDomain = _1;
template <typename ParticleType>
VTKM_EXEC void PostProcessStep(const vtkm::Id& vtkmNotUsed(idx),
ParticleType& vtkmNotUsed(particle)) const
{
}
template <typename ParticleType>
VTKM_EXEC void StepUpdate(const vtkm::Id& idx,
ParticleType& particle,
const vtkm::Vec3f& pos,
const vtkm::FloatDefault& time) const
{
particle.Pos = pos;
particle.NumSteps++;
particle.Time = time;
this->PostProcessStep(idx, particle);
}
template <typename ParticleType>
VTKM_EXEC bool ExtraCanContinue(const vtkm::Id& vtkmNotUsed(idx),
ParticleType& vtkmNotUsed(particle)) const
{
return true;
}
template <typename ParticleType>
VTKM_EXEC bool CanContinue(const vtkm::Id& idx, const ParticleType& particle) const
{
return (particle.Status.CheckOk() && !particle.Status.CheckTerminate() &&
!particle.Status.CheckSpatialBounds() && !particle.Status.CheckTemporalBounds() &&
this->ExtraCanContinue(idx, particle));
}
template <typename ParticleType>
VTKM_EXEC void UpdateStatus(const vtkm::Id& vtkmNotUsed(idx),
ParticleType& particle,
const IntegratorStatus& iStatus,
const vtkm::Id& maxSteps) const
{
if (particle.NumSteps == maxSteps)
particle.Status.SetTerminate();
if (iStatus.CheckFail())
particle.Status.SetFail();
if (iStatus.CheckSpatialBounds())
particle.Status.SetSpatialBounds();
if (iStatus.CheckTemporalBounds())
particle.Status.SetTemporalBounds();
}
template <typename IntegratorType, typename ParticleType>
VTKM_EXEC void operator()(const IntegratorType* integrator,
VTKM_EXEC void operator()(const vtkm::Id& idx,
const IntegratorType* integrator,
ParticleType& particle,
const vtkm::Id& maxSteps) const
@ -118,71 +176,47 @@ public:
vtkm::FloatDefault time = particle.Time;
IntegratorStatus status;
bool tookAnySteps = false;
while ((particle.Status & ParticleStatus::SUCCESS) != 0)
//the integrator status needs to be more robust:
// 1. you could have success AND at temporal boundary.
// 2. could you have success AND at spatial?
// 3. all three?
// std::cout<<idx<<": "<<inpos<<" #"<<particle.NumSteps<<" "<<particle.Status<<std::endl;
do
{
status = integrator->Step(inpos, time, outpos);
// If the status is OK, we only need to check if the particle
// has completed the maximum steps required.
if (status == IntegratorStatus::SUCCESS)
if (status.CheckOk())
{
//particle.TakeStep(idx, outpos, time);
particle.Pos = outpos;
particle.NumSteps++;
particle.Time = time;
if (particle.NumSteps == maxSteps)
{
particle.Status &= ~ParticleStatus::SUCCESS;
particle.Status |= ParticleStatus::TERMINATED;
}
// This is to keep track of the particle's time.
// This is what the Evaluator uses to determine if the particle
// has exited temporal boundary.
inpos = outpos;
this->StepUpdate(idx, particle, outpos, time);
tookAnySteps = true;
inpos = outpos;
}
// If the particle is at spatial or temporal boundary, take steps to just
// push it a little out of the boundary so that it will start advection in
// another domain, or in another time slice. Taking small steps enables
// reducing the error introduced at spatial or temporal boundaries.
else if (status == IntegratorStatus::OUTSIDE_TEMPORAL_BOUNDS)
{
particle.Status &= ~ParticleStatus::SUCCESS;
particle.Status |= ParticleStatus::EXIT_TEMPORAL_BOUNDARY;
break;
}
else if (status == IntegratorStatus::OUTSIDE_SPATIAL_BOUNDS)
{
status = integrator->SmallStep(inpos, time, outpos);
//particle.TakeStep(idx, outpos, time);
particle.Pos = outpos;
particle.NumSteps++;
particle.Time = time;
if (status == IntegratorStatus::OUTSIDE_TEMPORAL_BOUNDS)
//We can't take a step inside spatial boundary.
//Try and take a step just past the boundary.
else if (status.CheckSpatialBounds())
{
IntegratorStatus status2 = integrator->SmallStep(inpos, time, outpos);
if (status2.CheckOk())
{
particle.Status &= ~ParticleStatus::SUCCESS;
particle.Status |= ParticleStatus::EXIT_TEMPORAL_BOUNDARY;
break;
}
else if (status == IntegratorStatus::OUTSIDE_SPATIAL_BOUNDS)
{
particle.Status &= ~ParticleStatus::SUCCESS;
particle.Status |= ParticleStatus::EXIT_SPATIAL_BOUNDARY;
break;
}
else if (status == IntegratorStatus::FAIL)
{
particle.Status &= ~ParticleStatus::SUCCESS;
particle.Status |= ParticleStatus::FAIL;
break;
this->StepUpdate(idx, particle, outpos, time);
tookAnySteps = true;
//we took a step, so use this status to consider below.
status = status2;
}
}
}
this->UpdateStatus(idx, particle, status, maxSteps);
} while (this->CanContinue(idx, particle));
//Mark if any steps taken
if (tookAnySteps)
particle.Status |= ParticleStatus::TOOK_ANY_STEPS;
particle.Status.SetTookAnySteps();
else
particle.Status &= ~ParticleStatus::TOOK_ANY_STEPS;
particle.Status.ClearTookAnySteps();
}
};
@ -195,10 +229,11 @@ public:
using InputDomain = _1;
template <typename IntegratorType, typename IntegralCurveType>
VTKM_EXEC void operator()(const vtkm::Id& idx,
const IntegratorType* integrator,
IntegralCurveType& integralCurve) const
VTKM_EXEC void operator()(const vtkm::Id& /*idx*/,
const IntegratorType* /*integrator*/,
IntegralCurveType& /*integralCurve*/) const
{
#if 0
vtkm::Vec3f inpos = integralCurve.GetPos(idx);
vtkm::Vec3f outpos;
vtkm::FloatDefault time = integralCurve.GetTime(idx);
@ -251,6 +286,7 @@ public:
}
}
integralCurve.SetTookAnySteps(idx, tookAnySteps);
#endif
}
};
@ -293,34 +329,8 @@ public:
//AOS version
void Run(const IntegratorType& integrator,
vtkm::cont::ArrayHandle<vtkm::Particle>& p,
vtkm::cont::ArrayHandle<vtkm::Particle>& particles,
const vtkm::Id& MaxSteps)
{
using ParticleAdvectWorkletType = vtkm::worklet::particleadvection::ParticleAdvectWorkletAOS;
using ParticleWorkletDispatchType =
typename vtkm::worklet::DispatcherMapField<ParticleAdvectWorkletType>;
using ParticleTypeAOS = vtkm::worklet::particleadvection::ParticlesAOS;
vtkm::Id numSeeds = static_cast<vtkm::Id>(p.GetNumberOfValues());
//Create and invoke the particle advection.
vtkm::cont::ArrayHandleIndex idxArray(numSeeds);
ParticleTypeAOS particles(p, MaxSteps);
#ifdef VTKM_CUDA
// This worklet needs some extra space on CUDA.
vtkm::cont::cuda::ScopedCudaStackSize stack(16 * 1024);
(void)stack;
#endif // VTKM_CUDA
//Invoke particle advection worklet
ParticleWorkletDispatchType particleWorkletDispatch;
particleWorkletDispatch.Invoke(idxArray, integrator, particles);
}
//AOS version
void Run2(const IntegratorType& integrator,
vtkm::cont::ArrayHandle<vtkm::Particle>& particles,
const vtkm::Id& MaxSteps)
{
using ParticleAdvectWorkletType = vtkm::worklet::particleadvection::ParticleAdvectWorkletAOS2;
using ParticleWorkletDispatchType =
@ -329,6 +339,7 @@ public:
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::ArrayHandleIndex idxArray(numSeeds);
#ifdef VTKM_CUDA
// This worklet needs some extra space on CUDA.
@ -338,7 +349,7 @@ public:
//Invoke particle advection worklet
ParticleWorkletDispatchType particleWorkletDispatch;
particleWorkletDispatch.Invoke(integrator, particles, maxSteps);
particleWorkletDispatch.Invoke(idxArray, integrator, particles, maxSteps);
}
};

45
vtkm/worklet/particleadvection/Particles.h
View File

@ -70,66 +70,69 @@ public:
VTKM_EXEC
void SetOK(const vtkm::Id& idx)
{
Clear(idx);
Status.Set(idx, ParticleStatus::SUCCESS);
// Clear(idx);
// Status.Set(idx, ParticleStatus::SUCCESS);
}
VTKM_EXEC
void SetTerminated(const vtkm::Id& idx)
{
ClearBit(idx, ParticleStatus::SUCCESS);
SetBit(idx, ParticleStatus::TERMINATED);
// ClearBit(idx, ParticleStatus::SUCCESS);
// SetBit(idx, ParticleStatus::TERMINATED);
}
VTKM_EXEC
void SetTookAnySteps(const vtkm::Id& idx, const bool& val)
{
if (val)
SetBit(idx, ParticleStatus::TOOK_ANY_STEPS);
else
ClearBit(idx, ParticleStatus::TOOK_ANY_STEPS);
// if (val)
// SetBit(idx, ParticleStatus::TOOK_ANY_STEPS);
// else
// ClearBit(idx, ParticleStatus::TOOK_ANY_STEPS);
}
VTKM_EXEC
void SetExitSpatialBoundary(const vtkm::Id& idx)
{
ClearBit(idx, ParticleStatus::SUCCESS);
SetBit(idx, ParticleStatus::EXIT_SPATIAL_BOUNDARY);
// ClearBit(idx, ParticleStatus::SUCCESS);
// SetBit(idx, ParticleStatus::EXIT_SPATIAL_BOUNDARY);
}
VTKM_EXEC
void SetExitTemporalBoundary(const vtkm::Id& idx)
{
ClearBit(idx, ParticleStatus::SUCCESS);
SetBit(idx, ParticleStatus::EXIT_TEMPORAL_BOUNDARY);
// ClearBit(idx, ParticleStatus::SUCCESS);
// SetBit(idx, ParticleStatus::EXIT_TEMPORAL_BOUNDARY);
}
VTKM_EXEC
void SetError(const vtkm::Id& idx)
{
ClearBit(idx, ParticleStatus::SUCCESS);
SetBit(idx, ParticleStatus::FAIL);
// ClearBit(idx, ParticleStatus::SUCCESS);
// SetBit(idx, ParticleStatus::FAIL);
}
/* Check Status */
VTKM_EXEC
bool OK(const vtkm::Id& idx) { return CheckBit(idx, ParticleStatus::SUCCESS); }
bool OK(const vtkm::Id& idx) { return true; /*CheckBit(idx, ParticleStatus::SUCCESS);*/ }
VTKM_EXEC
bool Terminated(const vtkm::Id& idx) { return CheckBit(idx, ParticleStatus::TERMINATED); }
bool Terminated(const vtkm::Id& idx)
{
return true; /*CheckBit(idx, ParticleStatus::TERMINATED);*/
}
VTKM_EXEC
bool ExitSpatialBoundary(const vtkm::Id& idx)
{
return CheckBit(idx, ParticleStatus::EXIT_SPATIAL_BOUNDARY);
return true; //CheckBit(idx, ParticleStatus::EXIT_SPATIAL_BOUNDARY);
}
VTKM_EXEC
bool ExitTemporalBoundary(const vtkm::Id& idx)
{
return CheckBit(idx, ParticleStatus::EXIT_TEMPORAL_BOUNDARY);
return true; //CheckBit(idx, ParticleStatus::EXIT_TEMPORAL_BOUNDARY);
}
VTKM_EXEC
bool Error(const vtkm::Id& idx) { return CheckBit(idx, ParticleStatus::FAIL); }
bool Error(const vtkm::Id& idx) { return true; /*CheckBit(idx, ParticleStatus::FAIL);*/ }
VTKM_EXEC
bool Integrateable(const vtkm::Id& idx)
{
return OK(idx) && !(Terminated(idx) || ExitSpatialBoundary(idx) || ExitTemporalBoundary(idx));
return true; //OK(idx) && !(Terminated(idx) || ExitSpatialBoundary(idx) || ExitTemporalBoundary(idx));
}
VTKM_EXEC
bool Done(const vtkm::Id& idx) { return !Integrateable(idx); }
bool Done(const vtkm::Id& idx) { return true; /*!Integrateable(idx);*/ }
/* Bit Operations */
VTKM_EXEC