Commit 818782e3 authored by Abhishek Yenpure's avatar Abhishek Yenpure
Browse files

Adding a template for termination criteria

-- Adding a `Termination` template over particle to determine termination
-- Adding a NormalTermination which terminates a particles for duration as well as
   zero velocity
parent 3b8ead6c
......@@ -55,6 +55,10 @@ public:
VTKM_EXEC_CONT void ClearInGhostCell() { this->reset(this->IN_GHOST_CELL_BIT); }
VTKM_EXEC_CONT bool CheckInGhostCell() const { return this->test(this->IN_GHOST_CELL_BIT); }
VTKM_EXEC_CONT void SetZeroVelocity() { this->set(this->ZERO_VELOCITY); }
VTKM_EXEC_CONT void ClearZeroVelocity() { this->reset(this->ZERO_VELOCITY); }
VTKM_EXEC_CONT bool CheckZeroVelocity() const { return this->test(this->ZERO_VELOCITY); }
private:
static constexpr vtkm::Id SUCCESS_BIT = 0;
static constexpr vtkm::Id TERMINATE_BIT = 1;
......@@ -62,6 +66,7 @@ private:
static constexpr vtkm::Id TEMPORAL_BOUNDS_BIT = 3;
static constexpr vtkm::Id TOOK_ANY_STEPS_BIT = 4;
static constexpr vtkm::Id IN_GHOST_CELL_BIT = 5;
static constexpr vtkm::Id ZERO_VELOCITY = 6;
};
inline VTKM_CONT std::ostream& operator<<(std::ostream& s, const vtkm::ParticleStatus& status)
......@@ -71,6 +76,7 @@ inline VTKM_CONT std::ostream& operator<<(std::ostream& s, const vtkm::ParticleS
s << " spat= " << status.CheckSpatialBounds();
s << " temp= " << status.CheckTemporalBounds();
s << " ghst= " << status.CheckInGhostCell();
s << " zvol= " << status.CheckZeroVelocity();
s << "]";
return s;
}
......@@ -93,6 +99,8 @@ public:
, Status(status)
, Time(time)
{
this->Direction = 1;
this->DirChanged = false;
}
VTKM_EXEC_CONT
......@@ -100,6 +108,8 @@ public:
: Pos(p.Pos)
, ID(p.ID)
, NumSteps(p.NumSteps)
, Direction(p.Direction)
, DirChanged(p.DirChanged)
, Status(p.Status)
, Time(p.Time)
{
......@@ -140,6 +150,8 @@ public:
vtkm::Vec3f Pos;
vtkm::Id ID = -1;
vtkm::Id NumSteps = 0;
vtkm::Id Direction = 1;
bool DirChanged = 0;
vtkm::ParticleStatus Status;
vtkm::FloatDefault Time = 0;
};
......
......@@ -13,6 +13,8 @@
#ifndef vtk_m_worklet_particleadvection_EulerIntegrator_h
#define vtk_m_worklet_particleadvection_EulerIntegrator_h
#include <vtkm/VectorAnalysis.h>
namespace vtkm
{
namespace worklet
......@@ -38,10 +40,13 @@ public:
auto time = particle.Time;
auto inpos = particle.GetEvaluationPosition(stepLength);
vtkm::VecVariable<vtkm::Vec3f, 2> vectors;
GridEvaluatorStatus status = this->Evaluator.Evaluate(inpos, time, vectors);
if (status.CheckOk())
GridEvaluatorStatus evalStatus = this->Evaluator.Evaluate(inpos, time, vectors);
if (evalStatus.CheckOk())
velocity = particle.Velocity(vectors, stepLength);
return IntegratorStatus(status);
auto status = IntegratorStatus(evalStatus);
if(vtkm::MagnitudeSquared(velocity) <= vtkm::Epsilon<vtkm::FloatDefault>())
status.SetZeroVelocity();
return status;
}
private:
......
......@@ -38,12 +38,14 @@ public:
VTKM_EXEC_CONT IntegratorStatus(const bool& ok,
const bool& spatial,
const bool& temporal,
const bool& inGhost)
const bool& inGhost,
const bool& zero)
{
this->set(this->SUCCESS_BIT, ok);
this->set(this->SPATIAL_BOUNDS_BIT, spatial);
this->set(this->TEMPORAL_BOUNDS_BIT, temporal);
this->set(this->IN_GHOST_CELL_BIT, inGhost);
this->set(this->ZERO_VELOCITY, zero);
}
VTKM_EXEC_CONT IntegratorStatus(const GridEvaluatorStatus& es)
......@@ -69,17 +71,22 @@ public:
VTKM_EXEC_CONT void SetInGhostCell() { this->set(this->IN_GHOST_CELL_BIT); }
VTKM_EXEC_CONT bool CheckInGhostCell() const { return this->test(this->IN_GHOST_CELL_BIT); }
VTKM_EXEC_CONT void SetZeroVelocity() { this->set(this->ZERO_VELOCITY); }
VTKM_EXEC_CONT bool CheckZeroVelocity() const { return this->test(this->ZERO_VELOCITY); }
private:
static constexpr vtkm::Id SUCCESS_BIT = 0;
static constexpr vtkm::Id SPATIAL_BOUNDS_BIT = 1;
static constexpr vtkm::Id TEMPORAL_BOUNDS_BIT = 2;
static constexpr vtkm::Id IN_GHOST_CELL_BIT = 3;
static constexpr vtkm::Id ZERO_VELOCITY = 4;
};
inline VTKM_CONT std::ostream& operator<<(std::ostream& s, const IntegratorStatus& status)
{
s << "[ok= " << status.CheckOk() << " sp= " << status.CheckSpatialBounds()
<< " tm= " << status.CheckTemporalBounds() << " gc= " << status.CheckInGhostCell() << "]";
<< " tm= " << status.CheckTemporalBounds() << " gc= " << status.CheckInGhostCell()
<< " zv= " << status.CheckZeroVelocity() << "]";
return s;
}
}
......
......@@ -25,6 +25,7 @@
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/particleadvection/Particles.h>
#include <vtkm/worklet/particleadvection/Stepper.h>
#include <vtkm/worklet/particleadvection/Termination.h>
#ifdef VTKM_CUDA
#include <vtkm/cont/cuda/internal/ScopedCudaStackSize.h>
......@@ -53,6 +54,7 @@ public:
IntegralCurveType& integralCurve,
const vtkm::Id& maxSteps) const
{
(void)maxSteps;
auto particle = integralCurve.GetParticle(idx);
vtkm::FloatDefault time = particle.Time;
bool tookAnySteps = false;
......@@ -72,10 +74,9 @@ public:
integralCurve.StepUpdate(idx, particle, time, outpos);
tookAnySteps = true;
}
//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())
{
status = integrator.SmallStep(particle, time, outpos);
if (status.CheckOk())
......@@ -83,8 +84,8 @@ public:
integralCurve.StepUpdate(idx, particle, time, outpos);
tookAnySteps = true;
}
}
integralCurve.StatusUpdate(idx, status, maxSteps);
}*/
integralCurve.StatusUpdate(idx, status);
} while (integralCurve.CanContinue(idx));
//Mark if any steps taken
......@@ -109,7 +110,8 @@ public:
using ParticleAdvectWorkletType = vtkm::worklet::particleadvection::ParticleAdvectWorklet;
using ParticleWorkletDispatchType =
typename vtkm::worklet::DispatcherMapField<ParticleAdvectWorkletType>;
using ParticleArrayType = vtkm::worklet::particleadvection::Particles<ParticleType>;
using TerminationType = vtkm::worklet::particleadvection::NormalTermination<ParticleType>;
using ParticleArrayType = vtkm::worklet::particleadvection::Particles<ParticleType, TerminationType>;
vtkm::Id numSeeds = static_cast<vtkm::Id>(particles.GetNumberOfValues());
//Create and invoke the particle advection.
......@@ -126,8 +128,10 @@ public:
vtkm::cont::cuda::internal::ScopedCudaStackSize stack(16 * 1024);
(void)stack;
#endif // VTKM_CUDA
ParticleArrayType particlesObj(particles, MaxSteps);
TerminationType termination;
termination.SetMaxSteps(MaxSteps);
termination.SetZeroVelocityTermination(true);
ParticleArrayType particlesObj(particles, termination);
//Invoke particle advection worklet
ParticleWorkletDispatchType particleWorkletDispatch;
......@@ -188,8 +192,9 @@ public:
using ParticleWorkletDispatchType = typename vtkm::worklet::DispatcherMapField<
vtkm::worklet::particleadvection::ParticleAdvectWorklet>;
using TerminationType = vtkm::worklet::particleadvection::NormalTermination<ParticleType>;
using StreamlineArrayType =
vtkm::worklet::particleadvection::StateRecordingParticles<ParticleType>;
vtkm::worklet::particleadvection::StateRecordingParticles<ParticleType, TerminationType>;
vtkm::cont::ArrayHandle<vtkm::Id> initialStepsTaken;
......@@ -209,7 +214,10 @@ public:
#endif // VTKM_CUDA
//Run streamline worklet
StreamlineArrayType streamlines(particles, MaxSteps);
TerminationType termination;
termination.SetMaxSteps(MaxSteps);
termination.SetZeroVelocityTermination(true);
StreamlineArrayType streamlines(particles, MaxSteps, termination);
ParticleWorkletDispatchType particleWorkletDispatch;
vtkm::cont::ArrayHandleConstant<vtkm::Id> maxSteps(MaxSteps, numSeeds);
particleWorkletDispatch.Invoke(idxArray, it, streamlines, maxSteps);
......
......@@ -18,6 +18,7 @@
#include <vtkm/cont/ArrayHandleConstant.h>
#include <vtkm/cont/ExecutionObjectBase.h>
#include <vtkm/worklet/particleadvection/IntegratorStatus.h>
#include <vtkm/worklet/particleadvection/Termination.h>
namespace vtkm
{
......@@ -25,24 +26,24 @@ namespace worklet
{
namespace particleadvection
{
template <typename ParticleType>
template <typename ParticleType, typename TerminationType>
class ParticleExecutionObject
{
public:
VTKM_EXEC_CONT
ParticleExecutionObject()
: Particles()
, MaxSteps(0)
, Termination()
{
}
ParticleExecutionObject(vtkm::cont::ArrayHandle<ParticleType> particleArray,
vtkm::Id maxSteps,
const TerminationType& termination,
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token)
: Termination(termination)
{
Particles = particleArray.PrepareForInPlace(device, token);
MaxSteps = maxSteps;
}
VTKM_EXEC
......@@ -61,19 +62,16 @@ public:
newParticle.Pos = pt;
newParticle.Time = time;
newParticle.NumSteps++;
newParticle.DirChanged = false;
this->Particles.Set(idx, newParticle);
}
VTKM_EXEC
void StatusUpdate(const vtkm::Id& idx,
const vtkm::worklet::particleadvection::IntegratorStatus& status,
vtkm::Id maxSteps)
const vtkm::worklet::particleadvection::IntegratorStatus& status)
{
ParticleType p(this->GetParticle(idx));
if (p.NumSteps == maxSteps)
p.Status.SetTerminate();
if (status.CheckFail())
p.Status.SetFail();
if (status.CheckSpatialBounds())
......@@ -82,16 +80,19 @@ public:
p.Status.SetTemporalBounds();
if (status.CheckInGhostCell())
p.Status.SetInGhostCell();
if (status.CheckZeroVelocity())
p.Status.SetZeroVelocity();
this->Particles.Set(idx, p);
}
VTKM_EXEC
bool CanContinue(const vtkm::Id& idx)
{
ParticleType p(this->GetParticle(idx));
return (p.Status.CheckOk() && !p.Status.CheckTerminate() && !p.Status.CheckSpatialBounds() &&
!p.Status.CheckTemporalBounds() && !p.Status.CheckInGhostCell());
ParticleType particle(this->GetParticle(idx));
auto terminate = this->Termination.CheckTermination(particle);
this->Particles.Set(idx, particle);
return terminate;
}
VTKM_EXEC
......@@ -107,44 +108,46 @@ public:
protected:
using ParticlePortal = typename vtkm::cont::ArrayHandle<ParticleType>::WritePortalType;
ParticlePortal Particles;
vtkm::Id MaxSteps;
TerminationType Termination;
};
template <typename ParticleType>
template <typename ParticleType, typename TerminationType>
class Particles : public vtkm::cont::ExecutionObjectBase
{
protected:
vtkm::cont::ArrayHandle<ParticleType> ParticleArray;
TerminationType Termination;
public:
VTKM_CONT vtkm::worklet::particleadvection::ParticleExecutionObject<ParticleType>
PrepareForExecution(vtkm::cont::DeviceAdapterId device, vtkm::cont::Token& token) const
VTKM_CONT auto PrepareForExecution(vtkm::cont::DeviceAdapterId device, vtkm::cont::Token& token) const
-> vtkm::worklet::particleadvection::ParticleExecutionObject<ParticleType, decltype(this->Termination.PrepareForExecution(device, token))>
{
return vtkm::worklet::particleadvection::ParticleExecutionObject<ParticleType>(
this->ParticleArray, this->MaxSteps, device, token);
auto termination = this->Termination.PrepareForExecution(device, token);
return vtkm::worklet::particleadvection::ParticleExecutionObject<ParticleType, decltype(termination)>(
this->ParticleArray, termination, device, token);
}
VTKM_CONT
Particles(vtkm::cont::ArrayHandle<ParticleType>& pArray, vtkm::Id& maxSteps)
Particles(vtkm::cont::ArrayHandle<ParticleType>& pArray,
TerminationType& termination)
: ParticleArray(pArray)
, MaxSteps(maxSteps)
, Termination(termination)
{
}
Particles() {}
protected:
vtkm::cont::ArrayHandle<ParticleType> ParticleArray;
vtkm::Id MaxSteps;
};
template <typename ParticleType>
class StateRecordingParticleExecutionObject : public ParticleExecutionObject<ParticleType>
template <typename ParticleType, typename TerminationType>
class StateRecordingParticleExecutionObject : public ParticleExecutionObject<ParticleType, TerminationType>
{
public:
VTKM_EXEC_CONT
StateRecordingParticleExecutionObject()
: ParticleExecutionObject<ParticleType>()
: ParticleExecutionObject<ParticleType, TerminationType>()
, History()
, Length(0)
, StepCount()
......@@ -157,9 +160,10 @@ public:
vtkm::cont::ArrayHandle<vtkm::Id> validPointArray,
vtkm::cont::ArrayHandle<vtkm::Id> stepCountArray,
vtkm::Id maxSteps,
const TerminationType& termination,
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token)
: ParticleExecutionObject<ParticleType>(pArray, maxSteps, device, token)
: ParticleExecutionObject<ParticleType, TerminationType>(pArray, termination, device, token)
, Length(maxSteps + 1)
{
vtkm::Id numPos = pArray.GetNumberOfValues();
......@@ -171,7 +175,7 @@ public:
VTKM_EXEC
void PreStepUpdate(const vtkm::Id& idx)
{
ParticleType p = this->ParticleExecutionObject<ParticleType>::GetParticle(idx);
ParticleType p = this->ParticleExecutionObject<ParticleType, TerminationType>::GetParticle(idx);
if (this->StepCount.Get(idx) == 0)
{
vtkm::Id loc = idx * Length;
......@@ -187,7 +191,7 @@ public:
vtkm::FloatDefault time,
const vtkm::Vec3f& pt)
{
this->ParticleExecutionObject<ParticleType>::StepUpdate(idx, particle, time, pt);
this->ParticleExecutionObject<ParticleType, TerminationType>::StepUpdate(idx, particle, time, pt);
//local step count.
vtkm::Id stepCount = this->StepCount.Get(idx);
vtkm::Id loc = idx * Length + stepCount;
......@@ -206,9 +210,17 @@ protected:
IdPortal ValidPoint;
};
template <typename ParticleType>
template <typename ParticleType, typename TerminationType>
class StateRecordingParticles : vtkm::cont::ExecutionObjectBase
{
protected:
vtkm::cont::ArrayHandle<ParticleType> ParticleArray;
vtkm::cont::ArrayHandle<vtkm::Vec3f> HistoryArray;
vtkm::cont::ArrayHandle<vtkm::Id> StepCountArray;
vtkm::cont::ArrayHandle<vtkm::Id> ValidPointArray;
vtkm::Id MaxSteps;
TerminationType Termination;
public:
//Helper functor for compacting history
struct IsOne
......@@ -221,22 +233,26 @@ public:
};
VTKM_CONT vtkm::worklet::particleadvection::StateRecordingParticleExecutionObject<ParticleType>
PrepareForExecution(vtkm::cont::DeviceAdapterId device, vtkm::cont::Token& token) const
VTKM_CONT VTKM_CONT auto PrepareForExecution(vtkm::cont::DeviceAdapterId device, vtkm::cont::Token& token) const
-> vtkm::worklet::particleadvection::StateRecordingParticleExecutionObject<ParticleType, decltype(this->Termination.PrepareForExecution(device, token))>
{
return vtkm::worklet::particleadvection::StateRecordingParticleExecutionObject<ParticleType>(
auto termination = this->Termination.PrepareForExecution(device, token);
return vtkm::worklet::particleadvection::StateRecordingParticleExecutionObject<ParticleType, decltype(termination)>(
this->ParticleArray,
this->HistoryArray,
this->ValidPointArray,
this->StepCountArray,
this->MaxSteps,
termination,
device,
token);
}
VTKM_CONT
StateRecordingParticles(vtkm::cont::ArrayHandle<ParticleType>& pArray, const vtkm::Id& maxSteps)
: MaxSteps(maxSteps)
, ParticleArray(pArray)
StateRecordingParticles(vtkm::cont::ArrayHandle<ParticleType>& pArray, vtkm::Id& maxSteps, const TerminationType& termination)
: ParticleArray(pArray)
, MaxSteps(maxSteps)
, Termination(termination)
{
vtkm::Id numParticles = static_cast<vtkm::Id>(pArray.GetNumberOfValues());
......@@ -253,11 +269,13 @@ public:
StateRecordingParticles(vtkm::cont::ArrayHandle<ParticleType>& pArray,
vtkm::cont::ArrayHandle<vtkm::Vec3f>& historyArray,
vtkm::cont::ArrayHandle<vtkm::Id>& validPointArray,
const TerminationType& termination,
vtkm::Id& maxSteps)
{
ParticleArray = pArray;
HistoryArray = historyArray;
ValidPointArray = validPointArray;
Termination = termination;
MaxSteps = maxSteps;
}
......@@ -266,13 +284,6 @@ public:
{
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;
};
......
......@@ -13,6 +13,8 @@
#ifndef vtk_m_worklet_particleadvection_RK4Integrator_h
#define vtk_m_worklet_particleadvection_RK4Integrator_h
#include <vtkm/VectorAnalysis.h>
namespace vtkm
{
namespace worklet
......@@ -35,8 +37,11 @@ public:
vtkm::FloatDefault stepLength,
vtkm::Vec3f& velocity) const
{
auto time = particle.Time;
auto inpos = particle.GetEvaluationPosition(stepLength);
auto time = particle.Time;
auto inpos = particle.GetEvaluationPosition(stepLength);
auto dir = particle.Direction;
auto evalpos = inpos;
vtkm::FloatDefault boundary = this->Evaluator.GetTemporalBoundary(static_cast<vtkm::Id>(1));
if ((time + stepLength + vtkm::Epsilon<vtkm::FloatDefault>() - boundary) > 0.0)
stepLength = boundary - time;
......@@ -57,29 +62,35 @@ public:
GridEvaluatorStatus evalStatus;
evalStatus = this->Evaluator.Evaluate(inpos, time, k1);
evalStatus = this->Evaluator.Evaluate(evalpos, time, k1);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
v1 = particle.Velocity(k1, stepLength);
evalStatus = this->Evaluator.Evaluate(inpos + var1 * v1, var2, k2);
evalpos = inpos + (dir * var1 * v1);
evalStatus = this->Evaluator.Evaluate(evalpos, var2, k2);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
v2 = particle.Velocity(k2, stepLength);
evalStatus = this->Evaluator.Evaluate(inpos + var1 * v2, var2, k3);
evalpos = inpos + (dir * var1 * v2);
evalStatus = this->Evaluator.Evaluate(evalpos, var2, k3);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
v3 = particle.Velocity(k3, stepLength);
evalStatus = this->Evaluator.Evaluate(inpos + stepLength * v3, var3, k4);
evalpos = inpos + (dir * stepLength * v3);
evalStatus = this->Evaluator.Evaluate(evalpos, var3, k4);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
v4 = particle.Velocity(k4, stepLength);
velocity = (v1 + 2 * v2 + 2 * v3 + v4) / static_cast<vtkm::FloatDefault>(6);
return IntegratorStatus(evalStatus);
auto status = IntegratorStatus(evalStatus);
if(vtkm::MagnitudeSquared(velocity) <= vtkm::Epsilon<vtkm::FloatDefault>())
status.SetZeroVelocity();
return status;
}
private:
......
......@@ -64,7 +64,7 @@ public:
auto status = this->Integrator.CheckStep(particle, this->DeltaT, velocity);
if (status.CheckOk())
{
outpos = particle.Pos + this->DeltaT * velocity;
outpos = particle.Pos + (particle.Direction * this->DeltaT * velocity);
time += this->DeltaT;
}
else
......
#ifndef vtkm_worklet_particleadvection_termination
#define vtkm_worklet_particleadvection_termination
namespace vtkm
{
namespace worklet
{
namespace particleadvection
{
template <typename ParticleType>
class NormalTerminationExec
{
public:
VTKM_EXEC_CONT
NormalTerminationExec()
: MaxSteps(0)
, ZeroVelocityTerminate(0)
{
}
VTKM_EXEC_CONT
NormalTerminationExec(vtkm::Id maxSteps,
bool zeroVelocityTerminate)
: MaxSteps(maxSteps)
, ZeroVelocityTerminate(zeroVelocityTerminate)
{
}
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.Status;
if(particle.NumSteps == this->MaxSteps)
{
status.SetTerminate();
}
bool terminate = status.CheckOk()
&& !status.CheckTerminate()
&& !status.CheckSpatialBounds()
&& !status.CheckTemporalBounds()
&& !status.CheckInGhostCell();
if(this->ZeroVelocityTerminate)
{
terminate = terminate && !status.CheckZeroVelocity();
}
return terminate;
}
private:
vtkm::Id MaxSteps;
bool ZeroVelocityTerminate = 0;
};
template <typename ParticleType>
class NormalTermination : public vtkm::cont::ExecutionObjectBase
{
public:
VTKM_CONT
NormalTermination()
: MaxSteps(0)
, ZeroVelocityTerminate(0)
{
}
VTKM_CONT void SetMaxSteps(vtkm::Id steps){this->MaxSteps = steps;}