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

Add support for ghost cells in particle advection.

Browse Source
This commit is contained in:
Dave Pugmire 2020-10-19 16:46:48 -04:00
parent aa21bf4f94
commit 35627b287e
10 changed files with 231 additions and 163 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

73
vtkm/Particle.h
View File

@ -50,79 +50,30 @@ public:
VTKM_EXEC_CONT void ClearTookAnySteps() { this->reset(this->TOOK_ANY_STEPS_BIT); }
VTKM_EXEC_CONT bool CheckTookAnySteps() const { return this->test(this->TOOK_ANY_STEPS_BIT); }
VTKM_EXEC_CONT void SetInGhostCell() { this->set(this->IN_GHOST_CELL_BIT); }
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); }
private:
static constexpr vtkm::Id SUCCESS_BIT = 0;
static constexpr vtkm::Id TERMINATE_BIT = 1;
static constexpr vtkm::Id SPATIAL_BOUNDS_BIT = 2;
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;
};
inline VTKM_CONT std::ostream& operator<<(std::ostream& s, const vtkm::ParticleStatus& status)
{
s << "[" << status.CheckOk() << " " << status.CheckTerminate() << " "
<< status.CheckSpatialBounds() << " " << status.CheckTemporalBounds() << "]";
s << "[ok= " << status.CheckOk();
s << " term= " << status.CheckTerminate();
s << " spat= " << status.CheckSpatialBounds();
s << " temp= " << status.CheckTemporalBounds();
s << " ghst= " << status.CheckInGhostCell();
s << "]";
return s;
}
class ParticleBase
{
public:
VTKM_EXEC_CONT
ParticleBase() {}
VTKM_EXEC_CONT virtual ~ParticleBase() noexcept
{
// This must not be defaulted, since defaulted virtual destructors are
// troublesome with CUDA __host__ __device__ markup.
}
VTKM_EXEC_CONT
ParticleBase(const vtkm::Vec3f& p,
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)
, Status(status)
, Time(time)
{
}
VTKM_EXEC_CONT
ParticleBase(const vtkm::ParticleBase& p)
: Pos(p.Pos)
, ID(p.ID)
, NumSteps(p.NumSteps)
, Status(p.Status)
, Time(p.Time)
{
}
vtkm::ParticleBase& operator=(const vtkm::ParticleBase&) = default;
// The basic particle is only meant to be advected in a velocity
// field. In that case it is safe to assume that the velocity value
// will always be stored in the first location of vectors
VTKM_EXEC_CONT
virtual vtkm::Vec3f Next(const vtkm::VecVariable<vtkm::Vec3f, 2>&, const vtkm::FloatDefault&) = 0;
// The basic particle is only meant to be advected in a velocity
// field. In that case it is safe to assume that the velocity value
// will always be stored in the first location of vectors
VTKM_EXEC_CONT
virtual vtkm::Vec3f Velocity(const vtkm::VecVariable<vtkm::Vec3f, 2>&,
const vtkm::FloatDefault&) = 0;
vtkm::Vec3f Pos;
vtkm::Id ID = -1;
vtkm::Id NumSteps = 0;
vtkm::ParticleStatus Status;
vtkm::FloatDefault Time = 0;
};
class Particle
{
public:
@ -186,7 +137,7 @@ public:
vtkm::FloatDefault Time = 0;
};
class Electron : public vtkm::ParticleBase
class Electron
{
public:
VTKM_EXEC_CONT

3
vtkm/worklet/particleadvection/EulerIntegrator.h
View File

@ -67,7 +67,8 @@ public:
auto inpos = particle->Pos;
vtkm::VecVariable<vtkm::Vec3f, 2> vectors;
GridEvaluatorStatus status = this->Evaluator.Evaluate(inpos, time, vectors);
velocity = particle->Velocity(vectors, stepLength);
if (status.CheckOk())
velocity = particle->Velocity(vectors, stepLength);
return IntegratorStatus(status);
}
};

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

@ -37,12 +37,14 @@ class GridEvaluatorStatus : public vtkm::Bitset<vtkm::UInt8>
{
public:
VTKM_EXEC_CONT GridEvaluatorStatus(){};
VTKM_EXEC_CONT GridEvaluatorStatus(bool ok, bool spatial, bool temporal)
VTKM_EXEC_CONT GridEvaluatorStatus(bool ok, bool spatial, bool temporal, bool inGhost)
{
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);
};
VTKM_EXEC_CONT void SetOk() { this->set(this->SUCCESS_BIT); }
VTKM_EXEC_CONT bool CheckOk() const { return this->test(this->SUCCESS_BIT); }
@ -55,10 +57,14 @@ public:
VTKM_EXEC_CONT void SetTemporalBounds() { this->set(this->TEMPORAL_BOUNDS_BIT); }
VTKM_EXEC_CONT bool CheckTemporalBounds() const { return this->test(this->TEMPORAL_BOUNDS_BIT); }
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); }
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;
};
}
}

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

@ -12,6 +12,7 @@
#define vtk_m_worklet_particleadvection_GridEvaluators_h
#include <vtkm/Bitset.h>
#include <vtkm/CellClassification.h>
#include <vtkm/Types.h>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/cont/ArrayHandle.h>
@ -38,6 +39,8 @@ namespace particleadvection
template <typename DeviceAdapter, typename FieldType>
class ExecutionGridEvaluator
{
using GhostCellArrayType = vtkm::cont::ArrayHandle<vtkm::UInt8>;
public:
VTKM_CONT
ExecutionGridEvaluator() = default;
@ -47,12 +50,15 @@ public:
std::shared_ptr<vtkm::cont::CellInterpolationHelper> interpolationHelper,
const vtkm::Bounds& bounds,
const FieldType& field,
const GhostCellArrayType& ghostCells,
vtkm::cont::Token& token)
: Bounds(bounds)
, Field(field.PrepareForExecution(DeviceAdapter(), token))
, GhostCells(ghostCells.PrepareForInput(DeviceAdapter(), token))
, HaveGhostCells(ghostCells.GetNumberOfValues() > 0)
, InterpolationHelper(interpolationHelper->PrepareForExecution(DeviceAdapter(), token))
, Locator(locator->PrepareForExecution(DeviceAdapter(), token))
{
Locator = locator->PrepareForExecution(DeviceAdapter(), token);
InterpolationHelper = interpolationHelper->PrepareForExecution(DeviceAdapter(), token);
Field = field.PrepareForExecution(DeviceAdapter(), token);
}
template <typename Point>
@ -62,11 +68,15 @@ public:
Point parametric;
Locator->FindCell(point, cellId, parametric);
return cellId != -1;
if (cellId == -1)
return false;
else
return !this->InGhostCell(cellId);
}
VTKM_EXEC
bool IsWithinTemporalBoundary(const vtkm::FloatDefault vtkmNotUsed(time)) const { return true; }
bool IsWithinTemporalBoundary(const vtkm::FloatDefault& vtkmNotUsed(time)) const { return true; }
VTKM_EXEC
vtkm::Bounds GetSpatialBoundary() const { return this->Bounds; }
@ -80,49 +90,68 @@ public:
}
template <typename Point>
VTKM_EXEC GridEvaluatorStatus Evaluate(const Point& pos,
vtkm::FloatDefault vtkmNotUsed(time),
VTKM_EXEC GridEvaluatorStatus Evaluate(const Point& point,
const vtkm::FloatDefault& time,
vtkm::VecVariable<Point, 2>& out) const
{
return this->Evaluate(pos, out);
}
template <typename Point>
VTKM_EXEC GridEvaluatorStatus Evaluate(const Point& point, vtkm::VecVariable<Point, 2>& out) const
{
vtkm::Id cellId;
Point parametric;
GridEvaluatorStatus status;
status.SetOk();
if (!this->IsWithinTemporalBoundary(time))
{
status.SetFail();
status.SetTemporalBounds();
}
Locator->FindCell(point, cellId, parametric);
if (cellId == -1)
{
status.SetFail();
status.SetSpatialBounds();
return status;
}
else if (this->InGhostCell(cellId))
{
status.SetFail();
status.SetInGhostCell();
status.SetSpatialBounds();
}
vtkm::UInt8 cellShape;
vtkm::IdComponent nVerts;
vtkm::VecVariable<vtkm::Id, 8> ptIndices;
vtkm::VecVariable<vtkm::Vec3f, 8> fieldValues;
InterpolationHelper->GetCellInfo(cellId, cellShape, nVerts, ptIndices);
//If initial checks ok, then do the evaluation.
if (status.CheckOk())
{
vtkm::UInt8 cellShape;
vtkm::IdComponent nVerts;
vtkm::VecVariable<vtkm::Id, 8> ptIndices;
vtkm::VecVariable<vtkm::Vec3f, 8> fieldValues;
InterpolationHelper->GetCellInfo(cellId, cellShape, nVerts, ptIndices);
this->Field->GetValue(ptIndices, nVerts, parametric, cellShape, out);
//for (vtkm::IdComponent i = 0; i < nVerts; i++)
// fieldValues.Append(Field->GetValue(ptIndices[i]));
this->Field->GetValue(ptIndices, nVerts, parametric, cellShape, out);
status.SetOk();
}
//vtkm::exec::CellInterpolate(fieldValues, parametric, cellShape, out);
status.SetOk();
return status;
}
private:
const vtkm::exec::CellLocator* Locator;
const vtkm::exec::CellInterpolationHelper* InterpolationHelper;
const vtkm::worklet::particleadvection::ExecutionField* Field;
VTKM_EXEC bool InGhostCell(const vtkm::Id& cellId) const
{
if (this->HaveGhostCells && cellId != -1)
return GhostCells.Get(cellId) == vtkm::CellClassification::GHOST;
return false;
}
using GhostCellPortal = typename vtkm::cont::ArrayHandle<vtkm::UInt8>::template ExecutionTypes<
DeviceAdapter>::PortalConst;
vtkm::Bounds Bounds;
const vtkm::worklet::particleadvection::ExecutionField* Field;
const GhostCellPortal GhostCells;
bool HaveGhostCells;
const vtkm::exec::CellInterpolationHelper* InterpolationHelper;
const vtkm::exec::CellLocator* Locator;
};
template <typename FieldType>
@ -135,16 +164,55 @@ public:
vtkm::cont::ArrayHandleCartesianProduct<AxisHandle, AxisHandle, AxisHandle>;
using Structured2DType = vtkm::cont::CellSetStructured<2>;
using Structured3DType = vtkm::cont::CellSetStructured<3>;
using GhostCellArrayType = vtkm::cont::ArrayHandle<vtkm::UInt8>;
VTKM_CONT
GridEvaluator() = default;
VTKM_CONT
GridEvaluator(const vtkm::cont::DataSet& dataSet, const FieldType& field)
: Bounds(dataSet.GetCoordinateSystem().GetBounds())
, Field(field)
, GhostCellArray()
{
this->InitializeLocator(dataSet.GetCoordinateSystem(), dataSet.GetCellSet());
if (dataSet.HasCellField("vtkmGhostCells"))
{
auto arr = dataSet.GetCellField("vtkmGhostCells").GetData();
if (arr.IsType<GhostCellArrayType>())
this->GhostCellArray = arr.Cast<GhostCellArrayType>();
else
throw vtkm::cont::ErrorInternal("vtkmGhostCells not of type vtkm::UInt8");
}
}
VTKM_CONT
GridEvaluator(const vtkm::cont::CoordinateSystem& coordinates,
const vtkm::cont::DynamicCellSet& cellset,
const FieldType& field)
: Field(field)
, Bounds(coordinates.GetBounds())
: Bounds(coordinates.GetBounds())
, Field(field)
, GhostCellArray()
{
this->InitializeLocator(coordinates, cellset);
}
template <typename DeviceAdapter>
VTKM_CONT ExecutionGridEvaluator<DeviceAdapter, FieldType> PrepareForExecution(
DeviceAdapter,
vtkm::cont::Token& token) const
{
return ExecutionGridEvaluator<DeviceAdapter, FieldType>(this->Locator,
this->InterpolationHelper,
this->Bounds,
this->Field,
this->GhostCellArray,
token);
}
private:
VTKM_CONT void InitializeLocator(const vtkm::cont::CoordinateSystem& coordinates,
const vtkm::cont::DynamicCellSet& cellset)
{
if (cellset.IsSameType(Structured2DType()) || cellset.IsSameType(Structured3DType()))
{
@ -203,20 +271,11 @@ public:
throw vtkm::cont::ErrorInternal("Unsupported cellset type.");
}
template <typename DeviceAdapter>
VTKM_CONT ExecutionGridEvaluator<DeviceAdapter, FieldType> PrepareForExecution(
DeviceAdapter,
vtkm::cont::Token& token) const
{
return ExecutionGridEvaluator<DeviceAdapter, FieldType>(
this->Locator, this->InterpolationHelper, this->Bounds, this->Field, token);
}
private:
std::shared_ptr<vtkm::cont::CellLocator> Locator;
std::shared_ptr<vtkm::cont::CellInterpolationHelper> InterpolationHelper;
FieldType Field;
vtkm::Bounds Bounds;
FieldType Field;
GhostCellArrayType GhostCellArray;
std::shared_ptr<vtkm::cont::CellInterpolationHelper> InterpolationHelper;
std::shared_ptr<vtkm::cont::CellLocator> Locator;
};
} //namespace particleadvection

70
vtkm/worklet/particleadvection/IntegratorBase.h
View File

@ -115,24 +115,8 @@ protected:
vtkm::FloatDefault& time,
vtkm::Vec3f& outpos) const override
{
// If particle is out of either spatial or temporal boundary to begin with,
// then return the corresponding status.
IntegratorStatus status;
if (!this->Evaluator.IsWithinSpatialBoundary(particle->Pos))
{
status.SetFail();
status.SetSpatialBounds();
return status;
}
if (!this->Evaluator.IsWithinTemporalBoundary(particle->Time))
{
status.SetFail();
status.SetTemporalBounds();
return status;
}
vtkm::Vec3f velocity;
status = CheckStep(particle, this->StepLength, velocity);
auto status = this->CheckStep(particle, this->StepLength, velocity);
if (status.CheckOk())
{
outpos = particle->Pos + this->StepLength * velocity;
@ -149,17 +133,6 @@ protected:
vtkm::FloatDefault& time,
vtkm::Vec3f& outpos) const override
{
if (!this->Evaluator.IsWithinSpatialBoundary(particle->Pos))
{
outpos = particle->Pos;
return IntegratorStatus(false, true, false);
}
if (!this->Evaluator.IsWithinTemporalBoundary(particle->Time))
{
outpos = particle->Pos;
return IntegratorStatus(false, false, true);
}
//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.
@ -170,7 +143,7 @@ protected:
//The binary search will be between {0, this->StepLength}
vtkm::FloatDefault stepRange[2] = { 0, this->StepLength };
vtkm::Vec3f currPos(particle->Pos), currVelocity;
vtkm::Vec3f currPos(particle->Pos);
vtkm::VecVariable<vtkm::Vec3f, 2> currValue, tmp;
auto evalStatus = this->Evaluator.Evaluate(currPos, particle->Time, currValue);
if (evalStatus.CheckFail())
@ -182,44 +155,57 @@ protected:
{
//Try a step midway between stepRange[0] and stepRange[1]
div *= 2;
vtkm::FloatDefault stepCurr = stepRange[0] + (this->StepLength / div);
vtkm::FloatDefault currStep = stepRange[0] + (this->StepLength / div);
//See if we can step by stepCurr.
IntegratorStatus status = this->CheckStep(particle, stepCurr, currVelocity);
//See if we can step by currStep
vtkm::Vec3f currVelocity;
IntegratorStatus status = this->CheckStep(particle, currStep, currVelocity);
if (status.CheckOk()) //Integration step succedded.
{
//See if this point is in/out.
auto newPos = particle->Pos + stepCurr * currVelocity;
evalStatus = this->Evaluator.Evaluate(newPos, particle->Time + stepCurr, tmp);
auto newPos = particle->Pos + currStep * currVelocity;
evalStatus = this->Evaluator.Evaluate(newPos, particle->Time + currStep, tmp);
if (evalStatus.CheckOk())
{
//Point still in. Update currPos and set range to {stepCurr, stepRange[1]}
//Point still in. Update currPos and set range to {currStep, stepRange[1]}
currPos = newPos;
stepRange[0] = stepCurr;
stepRange[0] = currStep;
}
else
{
//The step succedded, but the next point is outside.
//Step too long. Set range to: {stepRange[0], stepCurr} and continue.
stepRange[1] = stepCurr;
//Step too long. Set range to: {stepRange[0], currStep} and continue.
stepRange[1] = currStep;
}
}
else
{
//Step too long. Set range to: {stepRange[0], stepCurr} and continue.
stepRange[1] = stepCurr;
stepRange[1] = currStep;
}
}
evalStatus = this->Evaluator.Evaluate(currPos, particle->Time + stepRange[0], currValue);
if (evalStatus.CheckFail())
//The eval at Time + stepRange[0] better be *inside*
VTKM_ASSERT(evalStatus.CheckOk() && !evalStatus.CheckSpatialBounds());
if (evalStatus.CheckFail() || evalStatus.CheckSpatialBounds())
return IntegratorStatus(evalStatus);
//Update the position and time.
outpos = currPos + stepRange[1] * particle->Velocity(currValue, stepRange[1]);
time += stepRange[1];
return IntegratorStatus(
true, true, !this->Evaluator.IsWithinTemporalBoundary(particle->Time));
//Get the evaluation status for the point that is *just* outside of the data.
evalStatus = this->Evaluator.Evaluate(outpos, time, currValue);
//The eval should fail, and the point should be outside.
VTKM_ASSERT(evalStatus.CheckFail() && evalStatus.CheckSpatialBounds());
IntegratorStatus status(evalStatus);
status.SetOk(); //status is ok.
return status;
}
VTKM_EXEC

15
vtkm/worklet/particleadvection/IntegratorStatus.h
View File

@ -33,11 +33,15 @@ class IntegratorStatus : public vtkm::Bitset<vtkm::UInt8>
public:
VTKM_EXEC_CONT IntegratorStatus() {}
VTKM_EXEC_CONT IntegratorStatus(const bool& ok, const bool& spatial, const bool& temporal)
VTKM_EXEC_CONT IntegratorStatus(const bool& ok,
const bool& spatial,
const bool& temporal,
const bool& inGhost)
{
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);
}
VTKM_EXEC_CONT IntegratorStatus(const GridEvaluatorStatus& es)
@ -45,6 +49,7 @@ public:
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 void SetOk() { this->set(this->SUCCESS_BIT); }
@ -59,16 +64,20 @@ public:
VTKM_EXEC_CONT void SetTemporalBounds() { this->set(this->TEMPORAL_BOUNDS_BIT); }
VTKM_EXEC_CONT bool CheckTemporalBounds() const { return this->test(this->TEMPORAL_BOUNDS_BIT); }
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); }
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;
};
inline VTKM_CONT std::ostream& operator<<(std::ostream& s, const IntegratorStatus& status)
{
s << "[" << status.CheckOk() << " " << status.CheckSpatialBounds() << " "
<< status.CheckTemporalBounds() << "]";
s << "[ok= " << status.CheckOk() << " sp= " << status.CheckSpatialBounds()
<< " tm= " << status.CheckTemporalBounds() << " gc= " << status.CheckInGhostCell() << "]";
return s;
}
}

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

@ -54,7 +54,6 @@ public:
{
auto particle = integralCurve.GetParticle(idx);
//vtkm::Vec3f currPos = particle.Pos;
vtkm::FloatDefault time = particle.Time;
bool tookAnySteps = false;
@ -79,21 +78,15 @@ public:
//Try and take a step just past the boundary.
else if (status.CheckSpatialBounds())
{
IntegratorStatus status2 = integrator->SmallStep(&particle, time, outpos);
if (status2.CheckOk())
status = integrator->SmallStep(&particle, time, outpos);
if (status.CheckOk())
{
integralCurve.StepUpdate(idx, time, outpos);
particle.Pos = outpos;
tookAnySteps = true;
//we took a step, so use this status to consider below.
status = status2;
}
else
status =
IntegratorStatus(true, status2.CheckSpatialBounds(), status2.CheckTemporalBounds());
}
integralCurve.StatusUpdate(idx, status, maxSteps);
} while (integralCurve.CanContinue(idx));
//Mark if any steps taken

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

@ -77,6 +77,8 @@ public:
p.Status.SetSpatialBounds();
if (status.CheckTemporalBounds())
p.Status.SetTemporalBounds();
if (status.CheckInGhostCell())
p.Status.SetInGhostCell();
this->Particles.Set(idx, p);
}
@ -86,7 +88,7 @@ public:
ParticleType p = this->GetParticle(idx);
return (p.Status.CheckOk() && !p.Status.CheckTerminate() && !p.Status.CheckSpatialBounds() &&
!p.Status.CheckTemporalBounds());
!p.Status.CheckTemporalBounds() && !p.Status.CheckInGhostCell());
}
VTKM_EXEC

4
vtkm/worklet/particleadvection/RK4Integrator.h
View File

@ -107,7 +107,9 @@ public:
v4 = particle->Velocity(k4, stepLength);
velocity = (v1 + 2 * v2 + 2 * v3 + v4) / static_cast<vtkm::FloatDefault>(6);
return IntegratorStatus(true, false, evalStatus.CheckTemporalBounds());
return IntegratorStatus(evalStatus);
// return IntegratorStatus(true, false, evalStatus.CheckTemporalBounds(), evalStatus.CheckInGhostCell());
}
};

61
vtkm/worklet/testing/UnitTestParticleAdvection.cxx
View File

@ -16,6 +16,7 @@
#include <vtkm/cont/DataSetBuilderRectilinear.h>
#include <vtkm/cont/DataSetBuilderUniform.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/GhostCellClassify.h>
#include <vtkm/io/VTKDataSetReader.h>
#include <vtkm/worklet/ParticleAdvection.h>
#include <vtkm/worklet/particleadvection/EulerIntegrator.h>
@ -301,7 +302,7 @@ public:
vtkm::Vec3f& pointOut) const
{
vtkm::VecVariable<vtkm::Vec3f, 2> values;
status = evaluator.Evaluate(pointIn.Pos, values);
status = evaluator.Evaluate(pointIn.Pos, pointIn.Time, values);
pointOut = values[0];
}
};
@ -411,6 +412,12 @@ void ValidateIntegratorForBoundary(const vtkm::Bounds& bounds,
for (vtkm::Id index = 0; index < numPoints; index++)
{
Status status = statusPortal.Get(index);
if (!status.CheckOk())
{
std::cout << __FILE__ << " " << __LINE__ << std::endl;
std::cout << __FILE__ << " " << __LINE__ << std::endl;
std::cout << "Fail coming... at " << index << std::endl;
}
VTKM_TEST_ASSERT(status.CheckOk(), "Error in evaluator for " + msg);
VTKM_TEST_ASSERT(status.CheckSpatialBounds(), "Error in evaluator for " + msg);
//Result should be push just outside of the bounds.
@ -516,6 +523,56 @@ void TestEvaluators()
}
}
void TestGhostCellEvaluators()
{
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
vtkm::Id numX = 6;
vtkm::Id numY = 6;
vtkm::Id numZ = 6;
vtkm::Id3 dims(numX + 1, numY + 1, numZ + 1);
vtkm::filter::GhostCellClassify addGhost;
auto ds = vtkm::cont::DataSetBuilderUniform::Create(dims);
auto dsWithGhost = addGhost.Execute(ds);
vtkm::cont::ArrayHandle<vtkm::Vec3f> vecField;
vtkm::Vec3f vec(1, 0, 0);
CreateConstantVectorField(dims[0] * dims[1] * dims[2], vec, vecField);
dsWithGhost.AddPointField("vec", vecField);
FieldType velocities(vecField);
GridEvalType gridEval(dsWithGhost, velocities);
vtkm::FloatDefault stepSize = 0.1;
RK4Type rk4(gridEval, stepSize);
vtkm::worklet::ParticleAdvection pa;
std::vector<vtkm::Particle> seeds;
seeds.push_back(vtkm::Particle(vtkm::Vec3f(.5, .5, .5), 0));
seeds.push_back(vtkm::Particle(vtkm::Vec3f(4, 3, 3), 1));
seeds.push_back(vtkm::Particle(vtkm::Vec3f(5.5, 5.5, 5.5), 1));
seeds.push_back(vtkm::Particle(vtkm::Vec3f(.5, 3, 3), 1));
auto seedArray = vtkm::cont::make_ArrayHandle(seeds, vtkm::CopyFlag::Off);
auto res = pa.Run(rk4, seedArray, 10000);
auto posPortal = res.Particles.ReadPortal();
vtkm::Id numSeeds = seedArray.GetNumberOfValues();
for (vtkm::Id i = 0; i < numSeeds; i++)
{
const auto& p = posPortal.Get(i);
VTKM_TEST_ASSERT(p.Status.CheckSpatialBounds(), "Particle did not leave the dataset.");
VTKM_TEST_ASSERT(p.Status.CheckInGhostCell(), "Particle did not end up in ghost cell.");
}
// for (int i = 0; i < seeds.size(); i++)
// std::cout<<posPortal.Get(i).Pos<<" "<<posPortal.Get(i).Time<<" "<<posPortal.Get(i).Status<<std::endl;
}
void ValidateParticleAdvectionResult(
const vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>& res,
vtkm::Id nSeeds,
@ -968,6 +1025,8 @@ void TestParticleAdvection()
{
TestIntegrators();
TestEvaluators();
TestGhostCellEvaluators();
TestParticleStatus();
TestWorkletsBasic();
TestParticleWorkletsWithDataSetTypes();