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

Terminate particles when it hits zero-velocity.

Browse Source
This commit is contained in:
Dave Pugmire 2022-08-05 16:06:51 -04:00
parent f97e6179d6
commit a5993a647f
6 changed files with 90 additions and 40 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
vtkm/Particle.h
View File

@ -55,6 +55,16 @@ 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); }
VTKM_EXEC_CONT bool CanContinue() const
{
return this->CheckOk() && !this->CheckTerminate() && !this->CheckSpatialBounds() &&
!this->CheckTemporalBounds() && !this->CheckInGhostCell() && !this->CheckZeroVelocity();
}
private:
static constexpr vtkm::Id SUCCESS_BIT = 0;
static constexpr vtkm::Id TERMINATE_BIT = 1;
@ -62,6 +72,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 +82,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 << " zvel= " << status.CheckZeroVelocity();
s << "]";
return s;
}

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

@ -613,37 +613,60 @@ void TestParticleStatus()
const vtkm::Id3 dims(5, 5, 5);
vtkm::Id nElements = dims[0] * dims[1] * dims[2];
FieldHandle fieldArray;
CreateConstantVectorField(nElements, vtkm::Vec3f(1, 0, 0), fieldArray);
auto dataSets = vtkm::worklet::testing::CreateAllDataSets(bounds, dims, false);
for (auto& ds : dataSets)
auto vecs = { vtkm::Vec3f(1, 0, 0), vtkm::Vec3f(1, 1, 1), vtkm::Vec3f(0, 0, 0) };
for (auto vec : vecs)
{
using FieldType = vtkm::worklet::flow::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
FieldHandle fieldArray;
CreateConstantVectorField(nElements, vec, fieldArray);
vtkm::Id maxSteps = 1000;
vtkm::FloatDefault stepSize = 0.01f;
auto dataSets = vtkm::worklet::testing::CreateAllDataSets(bounds, dims, false);
for (auto& ds : dataSets)
{
ds.AddPointField("vec", fieldArray);
FieldType velocities(fieldArray);
using FieldType = vtkm::worklet::flow::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
GridEvalType eval(ds, velocities);
Stepper rk4(eval, stepSize);
vtkm::Id maxSteps = 1000;
vtkm::FloatDefault stepSize = 0.01f;
vtkm::worklet::flow::ParticleAdvection pa;
std::vector<vtkm::Particle> pts;
pts.push_back(vtkm::Particle(vtkm::Vec3f(.5, .5, .5), 0));
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);
auto portal = seedsArray.ReadPortal();
FieldType velocities(fieldArray);
bool tookStep0 = portal.Get(0).Status.CheckTookAnySteps();
bool tookStep1 = portal.Get(1).Status.CheckTookAnySteps();
VTKM_TEST_ASSERT(tookStep0 == true, "Particle failed to take any steps");
VTKM_TEST_ASSERT(tookStep1 == false, "Particle took a step when it should not have.");
GridEvalType eval(ds, velocities);
Stepper rk4(eval, stepSize);
vtkm::worklet::flow::ParticleAdvection pa;
std::vector<vtkm::Particle> pts;
pts.push_back(vtkm::Particle(vtkm::Vec3f(.5, .5, .5), 0));
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);
auto portal = seedsArray.ReadPortal();
bool tookStep0 = portal.Get(0).Status.CheckTookAnySteps();
bool tookStep1 = portal.Get(1).Status.CheckTookAnySteps();
bool isZero0 = portal.Get(0).Status.CheckZeroVelocity();
bool isZero1 = portal.Get(1).Status.CheckZeroVelocity();
if (vtkm::Magnitude(vec) > 0)
{
VTKM_TEST_ASSERT(tookStep0 == true, "Particle failed to take any steps");
VTKM_TEST_ASSERT(tookStep1 == false, "Particle took a step when it should not have.");
VTKM_TEST_ASSERT(isZero0 == false, "Particle in zero velocity when it should not be.");
VTKM_TEST_ASSERT(isZero1 == false, "Particle in zero velocity when it should not be.");
}
else
{
VTKM_TEST_ASSERT(tookStep0 == true, "Particle failed to take any steps");
VTKM_TEST_ASSERT(tookStep1 == false, "Particle took a step when it should not have.");
VTKM_TEST_ASSERT(isZero0 == true, "Particle in zero velocity when it should not be.");
VTKM_TEST_ASSERT(isZero1 == false, "Particle in zero velocity when it should not be.");
VTKM_TEST_ASSERT(portal.Get(0).NumSteps == 1, "Particle should have taken only 1 step.");
}
}
}
}
@ -882,13 +905,17 @@ void TestParticleAdvectionFile(const std::string& fname,
void TestParticleAdvection()
{
/*
TestIntegrators();
TestEvaluators();
TestGhostCellEvaluators();
*/
TestParticleStatus();
/*
TestWorkletsBasic();
TestParticleWorkletsWithDataSetTypes();
*/
//Fusion test.
std::vector<vtkm::Vec3f> fusionPts, fusionEndPts;

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

@ -44,7 +44,7 @@ public:
GridEvaluatorStatus status = this->Evaluator.Evaluate(inpos, time, vectors);
if (status.CheckOk())
velocity = particle.Velocity(vectors, stepLength);
return IntegratorStatus(status);
return IntegratorStatus(status, velocity);
}
private:

26
vtkm/filter/flow/worklet/IntegratorStatus.h
View File

@ -38,20 +38,30 @@ public:
VTKM_EXEC_CONT IntegratorStatus(const bool& ok,
const bool& spatial,
const bool& temporal,
const bool& inGhost)
const bool& inGhost,
const bool& isZero)
{
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_BIT, isZero);
}
VTKM_EXEC_CONT IntegratorStatus(const GridEvaluatorStatus& es)
: IntegratorStatus(es.CheckOk(),
es.CheckSpatialBounds(),
es.CheckTemporalBounds(),
es.CheckInGhostCell(),
false)
{
this->set(this->SUCCESS_BIT, es.CheckOk());
this->set(this->SPATIAL_BOUNDS_BIT, es.CheckSpatialBounds());
this->set(this->TEMPORAL_BOUNDS_BIT, es.CheckTemporalBounds());
this->set(this->IN_GHOST_CELL_BIT, es.CheckInGhostCell());
}
VTKM_EXEC_CONT IntegratorStatus(const GridEvaluatorStatus& es, const vtkm::Vec3f& velocity)
: IntegratorStatus(es)
{
this->set(this->ZERO_VELOCITY_BIT,
vtkm::MagnitudeSquared(velocity) <= vtkm::Epsilon<vtkm::FloatDefault>());
}
VTKM_EXEC_CONT void SetOk() { this->set(this->SUCCESS_BIT); }
@ -68,18 +78,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_BIT); }
VTKM_EXEC_CONT bool CheckZeroVelocity() const { return this->test(this->ZERO_VELOCITY_BIT); }
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_BIT = 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()
<< "zero= " << status.CheckZeroVelocity() << " ]";
return s;
}

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

@ -78,17 +78,14 @@ 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());
}
bool CanContinue(const vtkm::Id& idx) { return this->GetParticle(idx).Status.CanContinue(); }
VTKM_EXEC
void UpdateTookSteps(const vtkm::Id& idx, bool val)

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

@ -82,7 +82,7 @@ public:
velocity = (v1 + 2 * v2 + 2 * v3 + v4) / static_cast<vtkm::FloatDefault>(6);
return IntegratorStatus(evalStatus);
return IntegratorStatus(evalStatus, velocity);
}
private: