Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-19 18:45:43 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Merge topic 'fix_advection_error'

Browse Source 
7ce3d72d Merge branch 'master' of https://gitlab.kitware.com/vtk/vtk-m into fix_advection_error
44254f07 Allocating arrays for worklet execution
1f3f10c0 Moving memory allocations/initializations to a single location
cfb7e7b4 Particle advection changing method inputs
51ae5a0c Additional changes for proper particle advection
fca79d21 Error fixes for particle advection filter

Acked-by: Kitware Robot <kwrobot@kitware.com>
Merge-request: !922
This commit is contained in:
Abhishek Yenpure 2017-09-15 02:50:58 +00:00 committed by Kitware Robot
commit 8db5c1988c
5 changed files with 62 additions and 32 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
vtkm/worklet/ParticleAdvection.h
View File

@ -80,8 +80,12 @@ public:
vtkm::cont::ArrayHandleConstant<vtkm::Id> init(0, numSeeds);
stepsTaken.Allocate(numSeeds);
vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>::Copy(init, stepsTaken);
worklet.Run(it, pts, nSteps, status, stepsTaken);
vtkm::cont::ArrayHandleConstant<vtkm::Id> statusOK(static_cast<vtkm::Id>(1), numSeeds);
status.Allocate(numSeeds);
vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>::Copy(statusOK, status);
worklet.Run(it, pts, nSteps, status, stepsTaken);
//Create output.
ParticleAdvectionResult<FieldType> res(pts, status, stepsTaken);
return res;
@ -94,7 +98,7 @@ public:
ParticleAdvectionResult<FieldType> Run(
const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
vtkm::cont::ArrayHandle<vtkm::Id>& stepsTaken,
const vtkm::cont::ArrayHandle<vtkm::Id>& stepsAlreadyTaken,
const vtkm::Id& nSteps,
const DeviceAdapter&)
{
@ -102,9 +106,18 @@ public:
FieldType,
DeviceAdapter>
worklet;
vtkm::cont::ArrayHandle<vtkm::Id> status;
worklet.Run(it, pts, nSteps, status, stepsTaken);
vtkm::cont::ArrayHandle<vtkm::Id> stepsTaken, status;
vtkm::Id numSeeds = static_cast<vtkm::Id>(pts.GetNumberOfValues());
//Allocate status and steps arrays.
stepsTaken.Allocate(numSeeds);
vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>::Copy(stepsAlreadyTaken, stepsTaken);
vtkm::cont::ArrayHandleConstant<vtkm::Id> statusOK(static_cast<vtkm::Id>(1), numSeeds);
status.Allocate(numSeeds);
vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>::Copy(statusOK, status);
worklet.Run(it, pts, nSteps, status, stepsTaken);
//Create output.
ParticleAdvectionResult<FieldType> res(pts, status, stepsTaken);
return res;

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

@ -57,7 +57,7 @@ public:
}
VTKM_EXEC_CONT
void GetBoundary(vtkm::Vec<FieldType, 3> dir, vtkm::Vec<FieldType, 3> dirBounds) const
void GetBoundary(vtkm::Vec<FieldType, 3>& dir, vtkm::Vec<FieldType, 3>& dirBounds) const
{
dirBounds[0] = static_cast<FieldType>(dir[0] > 0 ? bounds.X.Max : bounds.X.Min);
dirBounds[1] = static_cast<FieldType>(dir[1] > 0 ? bounds.Y.Max : bounds.Y.Min);
@ -101,7 +101,7 @@ public:
}
VTKM_EXEC_CONT
void GetBoundary(vtkm::Vec<FieldType, 3> dir, vtkm::Vec<FieldType, 3> dirBounds) const
void GetBoundary(vtkm::Vec<FieldType, 3>& dir, vtkm::Vec<FieldType, 3>& dirBounds) const
{
dirBounds[0] = static_cast<FieldType>(dir[0] > 0 ? bounds.X.Max : bounds.X.Min);
dirBounds[1] = static_cast<FieldType>(dir[1] > 0 ? bounds.Y.Max : bounds.Y.Min);
@ -210,7 +210,7 @@ public:
}
VTKM_EXEC_CONT
void GetBoundary(vtkm::Vec<FieldType, 3> dir, vtkm::Vec<FieldType, 3> dirBounds) const
void GetBoundary(vtkm::Vec<FieldType, 3>& dir, vtkm::Vec<FieldType, 3>& dirBounds) const
{
dirBounds[0] = static_cast<FieldType>(dir[0] > 0 ? bounds.X.Max : bounds.X.Min);
dirBounds[1] = static_cast<FieldType>(dir[1] > 0 ? bounds.Y.Max : bounds.Y.Min);
@ -378,7 +378,7 @@ public:
}
VTKM_EXEC_CONT
void GetBoundary(vtkm::Vec<FieldType, 3> dir, vtkm::Vec<FieldType, 3> dirBounds) const
void GetBoundary(vtkm::Vec<FieldType, 3>& dir, vtkm::Vec<FieldType, 3>& dirBounds) const
{
dirBounds[0] = static_cast<FieldType>(dir[0] > 0 ? bounds.X.Max : bounds.X.Min);
dirBounds[1] = static_cast<FieldType>(dir[1] > 0 ? bounds.Y.Max : bounds.Y.Min);

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

@ -79,46 +79,59 @@ public:
}
VTKM_EXEC
FieldType GetEscapeStepLength(const vtkm::Vec<FieldType, 3>& inpos,
FieldType stepLength,
vtkm::Vec<FieldType, 3>& velocity) const
ParticleStatus GetEscapeStepLength(const vtkm::Vec<FieldType, 3>& inpos,
FieldType& stepLength,
vtkm::Vec<FieldType, 3>& velocity) const
{
this->CheckStep(inpos, stepLength, velocity);
ParticleStatus status = this->CheckStep(inpos, stepLength, velocity);
if (status != ParticleStatus::STATUS_OK)
{
stepLength += this->Tolerance;
return status;
}
FieldType magnitude = vtkm::Magnitude(velocity);
vtkm::Vec<FieldType, 3> dir = velocity / magnitude;
vtkm::Vec<FieldType, 3> dirBounds;
this->Evaluator.GetBoundary(dir, dirBounds);
/*Add a fraction just push the particle beyond the bounds*/
FieldType hx = (std::abs(dirBounds[0] - inpos[0]) + this->Tolerance) / std::abs(velocity[0]);
FieldType hy = (std::abs(dirBounds[1] - inpos[1]) + this->Tolerance) / std::abs(velocity[1]);
FieldType hz = (std::abs(dirBounds[2] - inpos[2]) + this->Tolerance) / std::abs(velocity[2]);
return vtkm::Min(hx, vtkm::Min(hy, hz));
FieldType hx = (vtkm::Abs(dirBounds[0] - inpos[0]) + this->Tolerance) / vtkm::Abs(velocity[0]);
FieldType hy = (vtkm::Abs(dirBounds[1] - inpos[1]) + this->Tolerance) / vtkm::Abs(velocity[1]);
FieldType hz = (vtkm::Abs(dirBounds[2] - inpos[2]) + this->Tolerance) / vtkm::Abs(velocity[2]);
stepLength = vtkm::Min(hx, vtkm::Min(hy, hz));
return status;
}
VTKM_EXEC
ParticleStatus PushOutOfDomain(vtkm::Vec<FieldType, 3> inpos,
ParticleStatus PushOutOfDomain(const vtkm::Vec<FieldType, 3>& inpos,
vtkm::Id numSteps,
vtkm::Vec<FieldType, 3>& outpos) const
{
ParticleStatus status;
outpos = inpos;
numSteps = (numSteps == 0) ? 1 : numSteps;
FieldType totalTime = static_cast<FieldType>(numSteps) * this->StepLength;
FieldType timeFraction = totalTime * this->Tolerance;
FieldType stepLength = this->StepLength / 2;
vtkm::Vec<FieldType, 3> velocity;
vtkm::Vec<FieldType, 3> velocity, currentVelocity;
this->CheckStep(inpos, 0.0f, currentVelocity);
if (this->ShortStepsSupported)
{
do
{
ParticleStatus status = this->CheckStep(inpos, stepLength, velocity);
status = this->CheckStep(inpos, stepLength, velocity);
if (status == ParticleStatus::STATUS_OK)
{
inpos = inpos + stepLength * velocity;
outpos = outpos + stepLength * velocity;
currentVelocity = velocity;
}
stepLength = stepLength / 2;
} while (stepLength > timeFraction);
}
stepLength = GetEscapeStepLength(inpos, stepLength, velocity);
outpos = inpos + stepLength * velocity;
status = GetEscapeStepLength(inpos, stepLength, velocity);
if (status != ParticleStatus::STATUS_OK)
outpos = outpos + stepLength * currentVelocity;
else
outpos = outpos + stepLength * velocity;
return ParticleStatus::EXITED_SPATIAL_BOUNDARY;
}

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

@ -114,20 +114,16 @@ private:
void run(vtkm::cont::ArrayHandle<vtkm::Id, FieldStorage>& statusArray,
vtkm::cont::ArrayHandle<vtkm::Id, FieldStorage>& stepsTaken)
{
typedef typename vtkm::worklet::DispatcherMapField<ParticleAdvectWorkletType>
ParticleWorkletDispatchType;
typedef vtkm::worklet::particleadvection::Particles<FieldType, DeviceAdapterTag> ParticleType;
using ParticleWorkletDispatchType =
typename vtkm::worklet::DispatcherMapField<ParticleAdvectWorkletType>;
using ParticleType = vtkm::worklet::particleadvection::Particles<FieldType, DeviceAdapterTag>;
vtkm::Id numSeeds = static_cast<vtkm::Id>(seedArray.GetNumberOfValues());
//Allocate status and steps arrays.
vtkm::cont::ArrayHandleConstant<vtkm::Id> ok(ParticleStatus::STATUS_OK, numSeeds);
statusArray.Allocate(numSeeds);
DeviceAlgorithm::Copy(ok, statusArray);
//Create and invoke the particle advection.
vtkm::cont::ArrayHandleIndex idxArray(numSeeds);
ParticleType particles(seedArray, stepsTaken, statusArray, maxSteps);
//Invoke particle advection worklet
ParticleAdvectWorkletType particleWorklet(integrator);
ParticleWorkletDispatchType particleWorkletDispatch(particleWorklet);
particleWorkletDispatch.Invoke(idxArray, particles);

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

@ -130,9 +130,17 @@ public:
SetBit(idx, TERMINATED);
}
VTKM_EXEC
void SetExitedSpatialBoundary(const vtkm::Id& idx) { SetBit(idx, EXITED_SPATIAL_BOUNDARY); }
void SetExitedSpatialBoundary(const vtkm::Id& idx)
{
ClearBit(idx, STATUS_OK);
SetBit(idx, EXITED_SPATIAL_BOUNDARY);
}
VTKM_EXEC
void SetExitedTemporalBoundary(const vtkm::Id& idx) { SetBit(idx, EXITED_TEMPORAL_BOUNDARY); }
void SetExitedTemporalBoundary(const vtkm::Id& idx)
{
ClearBit(idx, STATUS_OK);
SetBit(idx, EXITED_TEMPORAL_BOUNDARY);
}
VTKM_EXEC
void SetError(const vtkm::Id& idx)
{