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

Adding working temporal interpolation

Browse Source
This commit is contained in:
ayenpure 2018-03-23 00:49:32 -07:00
parent 8ae9c733e8
commit 077988aeb4
7 changed files with 132 additions and 179 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

164
examples/temporalparticleadvection/ParticleAdvection.cxx
View File

@ -34,40 +34,26 @@
#include <cstdlib>
#include <vector>
const vtkm::Id SPARSE = 0;
const vtkm::Id DENSE = 1;
const vtkm::Id MEDIUM = 2;
template <typename T>
static vtkm::Range subRange(vtkm::Range& range, T a, T b)
{
vtkm::Float32 arg1, arg2, len;
arg1 = static_cast<vtkm::Float32>(a);
arg2 = static_cast<vtkm::Float32>(b);
len = static_cast<vtkm::Float32>(range.Length());
return vtkm::Range(range.Min + arg1 * len, range.Min + arg2 * len);
}
template <typename T>
void ignore(T&&)
{
}
void RunTest(const std::string& fname,
vtkm::Id numSeeds,
vtkm::Id numSteps,
vtkm::Float32 stepSize,
vtkm::Id seeding)
void RunTest(vtkm::Id numSteps, vtkm::Float32 stepSize, vtkm::Id advectType)
{
using DeviceAdapter = VTKM_DEFAULT_DEVICE_ADAPTER_TAG;
using FieldType = vtkm::Float32;
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>>;
using FieldPortalConstType =
typename FieldHandle::template ExecutionTypes<DeviceAdapter>::PortalConst;
vtkm::io::reader::BOVDataSetReader reader(fname);
vtkm::cont::DataSet ds = reader.ReadDataSet();
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> fieldArray1;
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> fieldArray2;
vtkm::Id numValues;
vtkm::io::reader::BOVDataSetReader reader1("slice1.bov");
vtkm::cont::DataSet ds1 = reader1.ReadDataSet();
ds1.GetField(0).GetData().CopyTo(fieldArray1);
vtkm::io::reader::BOVDataSetReader reader2("slice2.bov");
vtkm::cont::DataSet ds2 = reader2.ReadDataSet();
ds2.GetField(0).GetData().CopyTo(fieldArray2);
using GridEvaluator =
vtkm::worklet::particleadvection::TemporalGridEvaluator<FieldPortalConstType,
@ -75,125 +61,67 @@ void RunTest(const std::string& fname,
DeviceAdapter>;
using Integrator = vtkm::worklet::particleadvection::EulerIntegrator<GridEvaluator, FieldType>;
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> fieldArray;
ds.GetField(0).GetData().CopyTo(fieldArray);
//GridEvaluator eval(ds.GetCoordinateSystem(), ds.GetCellSet(0), fieldArray);
GridEvaluator eval(ds.GetCoordinateSystem(),
ds.GetCellSet(0),
fieldArray,
GridEvaluator eval(ds1.GetCoordinateSystem(),
ds1.GetCellSet(0),
fieldArray1,
0,
ds.GetCoordinateSystem(),
ds.GetCellSet(0),
fieldArray,
1.0);
ds2.GetCoordinateSystem(),
ds2.GetCellSet(0),
fieldArray2,
10.0);
Integrator integrator(eval, stepSize);
std::vector<vtkm::Vec<FieldType, 3>> seeds;
srand(314);
vtkm::Bounds bounds = ds.GetCoordinateSystem().GetBounds();
if (seeding == SPARSE)
bounds = ds.GetCoordinateSystem().GetBounds();
else if (seeding == DENSE)
{
if (fname.find("astro") != std::string::npos)
{
bounds.X = subRange(bounds.X, .1, .15);
bounds.Y = subRange(bounds.Y, .1, .15);
bounds.Z = subRange(bounds.Z, .1, .15);
}
else if (fname.find("fusion") != std::string::npos)
{
bounds.X = subRange(bounds.X, .8, .85);
bounds.Y = subRange(bounds.Y, .55, .60);
bounds.Z = subRange(bounds.Z, .55, .60);
}
else if (fname.find("fishtank") != std::string::npos)
{
bounds.X = subRange(bounds.X, .1, .15);
bounds.Y = subRange(bounds.Y, .1, .15);
bounds.Z = subRange(bounds.Z, .55, .60);
}
}
else if (seeding == MEDIUM)
{
if (fname.find("astro") != std::string::npos)
{
bounds.X = subRange(bounds.X, .4, .6);
bounds.Y = subRange(bounds.Y, .4, .6);
bounds.Z = subRange(bounds.Z, .4, .6);
}
else if (fname.find("fusion") != std::string::npos)
{
bounds.X = subRange(bounds.X, .01, .99);
bounds.Y = subRange(bounds.Y, .01, .99);
bounds.Z = subRange(bounds.Z, .45, .55);
}
else if (fname.find("fishtank") != std::string::npos)
{
bounds.X = subRange(bounds.X, .4, .6);
bounds.Y = subRange(bounds.Y, .4, .6);
bounds.Z = subRange(bounds.Z, .4, .6);
}
}
for (int i = 0; i < numSeeds; i++)
vtkm::Id x = 0, y = 5, z = 0;
for (int i = 0; i <= 11; i++)
{
vtkm::Vec<FieldType, 3> point;
vtkm::Float32 rx = (vtkm::Float32)rand() / (vtkm::Float32)RAND_MAX;
vtkm::Float32 ry = (vtkm::Float32)rand() / (vtkm::Float32)RAND_MAX;
vtkm::Float32 rz = (vtkm::Float32)rand() / (vtkm::Float32)RAND_MAX;
point[0] = static_cast<FieldType>(bounds.X.Min + rx * bounds.X.Length());
point[1] = static_cast<FieldType>(bounds.Y.Min + ry * bounds.Y.Length());
point[2] = static_cast<FieldType>(bounds.Z.Min + rz * bounds.Z.Length());
point[0] = x;
point[1] = y;
point[2] = z++;
seeds.push_back(point);
}
/*#ifdef __BUILDING_TBB_VERSION__
int nT = tbb::task_scheduler_init::default_num_threads();
if (numThreads != -1)
nT = (int)numThreads;
//make sure the task_scheduler_init object is in scope when running sth w/ TBB
tbb::task_scheduler_init init(nT);
#else
ignore(numThreads);
#endif
*/
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> seedArray;
seedArray = vtkm::cont::make_ArrayHandle(seeds);
// if (advectType == 0)
// {
vtkm::worklet::ParticleAdvection particleAdvection;
particleAdvection.Run(integrator, seedArray, numSteps, DeviceAdapter());
// }
/*else
if (advectType == 0)
{
vtkm::worklet::ParticleAdvection particleAdvection;
particleAdvection.Run(integrator, seedArray, numSteps, DeviceAdapter());
}
else
{
vtkm::worklet::Streamline streamline;
streamline.Run(integrator, seedArray, numSteps, DeviceAdapter());
}*/
vtkm::worklet::StreamlineResult<FieldType> res =
streamline.Run(integrator, seedArray, numSteps, DeviceAdapter());
vtkm::cont::DataSet outData;
vtkm::cont::CoordinateSystem outputCoords("coordinates", res.positions);
outData.AddCellSet(res.polyLines);
outData.AddCoordinateSystem(outputCoords);
}
}
int main(int argc, char** argv)
{
vtkm::Id numSeeds, numSteps;
vtkm::Id numSteps;
vtkm::Float32 stepSize;
std::string dataFile;
vtkm::Id seeding = SPARSE;
vtkm::Id advectionType;
if (argc < 5)
if (argc < 4)
{
std::cout << "Wrong number of parameters provided" << std::endl;
exit(EXIT_FAILURE);
}
dataFile = std::string(argv[1]);
numSeeds = atoi(argv[2]);
numSteps = atoi(argv[3]);
stepSize = atof(argv[4]);
numSteps = atoi(argv[1]);
stepSize = atof(argv[2]);
advectionType = atoi(argv[3]);
RunTest(numSteps, stepSize, advectionType);
RunTest(dataFile, numSeeds, numSteps, stepSize, seeding);
return 0;
}

76
vtkm/worklet/ParticleAdvection.h
View File

@ -38,6 +38,7 @@ struct ParticleAdvectionResult
: positions()
, status()
, stepsTaken()
, times()
{
}
@ -50,9 +51,22 @@ struct ParticleAdvectionResult
{
}
ParticleAdvectionResult(const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>>& pos,
const vtkm::cont::ArrayHandle<vtkm::Id>& stat,
const vtkm::cont::ArrayHandle<vtkm::Id>& steps,
const vtkm::cont::ArrayHandle<FieldType>& timeArray)
: positions(pos)
, status(stat)
, stepsTaken(steps)
, times(timeArray)
{
}
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> positions;
vtkm::cont::ArrayHandle<vtkm::Id> status;
vtkm::cont::ArrayHandle<vtkm::Id> stepsTaken;
vtkm::cont::ArrayHandle<FieldType> times;
};
class ParticleAdvection
@ -76,7 +90,10 @@ public:
worklet;
vtkm::cont::ArrayHandle<vtkm::Id> stepsTaken, status;
vtkm::cont::ArrayHandle<FieldType> timeArray;
vtkm::Id numSeeds = static_cast<vtkm::Id>(pts.GetNumberOfValues());
//Allocate status and steps arrays.
vtkm::cont::ArrayHandleConstant<vtkm::Id> init(0, numSeeds);
stepsTaken.Allocate(numSeeds);
@ -86,41 +103,14 @@ public:
status.Allocate(numSeeds);
vtkm::cont::ArrayCopy(statusOK, status, DeviceAdapter());
worklet.Run(it, pts, nSteps, status, stepsTaken);
//Allocate memory to store the time for temporal integration.
vtkm::cont::ArrayHandleConstant<FieldType> time(0, numSeeds);
timeArray.Allocate(numSeeds);
vtkm::cont::ArrayCopy(time, timeArray, DeviceAdapter());
worklet.Run(it, pts, nSteps, status, stepsTaken, timeArray);
//Create output.
ParticleAdvectionResult<FieldType> res(pts, status, stepsTaken);
return res;
}
template <typename IntegratorType,
typename FieldType,
typename PointStorage,
typename DeviceAdapter>
ParticleAdvectionResult<FieldType> Run(
const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
const vtkm::cont::ArrayHandle<vtkm::Id>& stepsAlreadyTaken,
const vtkm::Id& nSteps,
const DeviceAdapter&)
{
vtkm::worklet::particleadvection::ParticleAdvectionWorklet<IntegratorType,
FieldType,
DeviceAdapter>
worklet;
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::ArrayCopy(stepsAlreadyTaken, stepsTaken, DeviceAdapter());
vtkm::cont::ArrayHandleConstant<vtkm::Id> statusOK(static_cast<vtkm::Id>(1), numSeeds);
status.Allocate(numSeeds);
vtkm::cont::ArrayCopy(statusOK, status, DeviceAdapter());
worklet.Run(it, pts, nSteps, status, stepsTaken);
//Create output.
ParticleAdvectionResult<FieldType> res(pts, status, stepsTaken);
ParticleAdvectionResult<FieldType> res(pts, status, stepsTaken, timeArray);
return res;
}
};
@ -177,18 +167,24 @@ public:
//Allocate status and steps arrays.
vtkm::Id numSeeds = seedArray.GetNumberOfValues();
vtkm::Id val = vtkm::worklet::particleadvection::ParticleStatus::STATUS_OK;
vtkm::cont::ArrayHandle<vtkm::Id> status, steps;
vtkm::cont::ArrayHandleConstant<vtkm::Id> ok(val, numSeeds);
status.Allocate(numSeeds);
DeviceAlgorithm::Copy(ok, status);
vtkm::cont::ArrayHandle<vtkm::Id> status, steps;
vtkm::cont::ArrayHandle<FieldType> timeArray;
vtkm::cont::ArrayHandleConstant<vtkm::Id> statusOK(static_cast<vtkm::Id>(1), numSeeds);
status.Allocate(numSeeds);
DeviceAlgorithm::Copy(statusOK, status);
vtkm::cont::ArrayHandleConstant<vtkm::Id> zero(0, numSeeds);
steps.Allocate(numSeeds);
DeviceAlgorithm::Copy(zero, steps);
worklet.Run(it, seedArray, nSteps, positions, polyLines, status, steps);
//Allocate memory to store the time for temporal integration.
vtkm::cont::ArrayHandleConstant<FieldType> time(0, numSeeds);
timeArray.Allocate(numSeeds);
DeviceAlgorithm::Copy(time, timeArray);
worklet.Run(it, seedArray, nSteps, positions, polyLines, status, steps, timeArray);
StreamlineResult<FieldType> res(positions, polyLines, status, steps);
return res;

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

@ -92,6 +92,9 @@ public:
}
}
VTKM_EXEC
FieldType GetStepLength() const { return StepLength; }
private:
FieldEvaluateType Evaluator;
FieldType StepLength;
@ -150,6 +153,9 @@ public:
return ParticleStatus::EXITED_SPATIAL_BOUNDARY;
}
VTKM_EXEC
FieldType GetStepLength() const { return StepLength; }
private:
FieldEvaluateType Evaluator;
FieldType StepLength;

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

@ -49,29 +49,37 @@ public:
template <typename IntegralCurveType>
VTKM_EXEC void operator()(const vtkm::Id& idx, IntegralCurveType& ic) const
{
vtkm::Vec<FieldType, 3> inpos = ic.GetPos(idx);
FieldType time = ic.GetTime(idx);
vtkm::Vec<FieldType, 3> outpos;
std::cout << "Advecting particle " << idx << std::endl;
vtkm::Vec<FieldType, 3> inpos = ic.GetPos(idx);
FieldType stepLength = integrator.GetStepLength();
vtkm::Vec<FieldType, 3> outpos;
std::cout << "Starting at : " << inpos[0] << ", " << inpos[1] << ", " << inpos[2] << std::endl;
while (!ic.Done(idx))
{
FieldType time = ic.GetTime(idx);
ParticleStatus status = integrator.Step(inpos, time, outpos);
if (status == ParticleStatus::STATUS_OK)
{
ic.TakeStep(idx, outpos, status);
time += stepLength;
std::cout << "Time : " << time << " at ";
ic.SetTime(idx, time);
std::cout << outpos[0] << ", " << outpos[1] << ", " << outpos[2] << std::endl;
inpos = outpos;
}
else if (status == ParticleStatus::EXITED_SPATIAL_BOUNDARY)
{
ic.TakeStep(idx, outpos, status);
//ic.TakeStep(idx, outpos, status);
ic.SetExitedSpatialBoundary(idx);
}
else if (status == ParticleStatus::EXITED_TEMPORAL_BOUNDARY)
{
ic.TakeStep(idx, outpos, status);
//ic.TakeStep(idx, outpos, status);
ic.SetExitedTemporalBoundary(idx);
}
}
std::cout << "Finished particle " << idx << std::endl;
}
ParticleAdvectWorklet(const IntegratorType& it)
@ -200,7 +208,7 @@ private:
vtkm::cont::ArrayHandleIndex idxArray(numSeeds);
vtkm::cont::ArrayHandle<vtkm::Id> validPoint;
std::vector<vtkm::Id> vpa(static_cast<std::size_t>(numSeeds * maxSteps), 0);
std::vector<vtkm::Id> vpa(numSeeds * maxSteps, 0);
validPoint = vtkm::cont::make_ArrayHandle(vpa);
//Compact history into positions.

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

@ -83,20 +83,31 @@ public:
VTKM_EXEC
void TakeStep(const vtkm::Id& idx, const vtkm::Vec<T, 3>& pt, ParticleStatus status)
{
// Irrespective of what the advected status of the particle is,
// we need to set the output position as the last step taken by
// the particle, and increase the number of steps take by 1.
Pos.Set(idx, pt);
vtkm::Id nSteps = Steps.Get(idx);
Steps.Set(idx, ++nSteps);
// If the status is OK, we only need to check if the particle
// has completed the maximum steps required.
if (status == ParticleStatus::STATUS_OK)
{
Pos.Set(idx, pt);
vtkm::Id nSteps = Steps.Get(idx);
nSteps = nSteps + 1;
Steps.Set(idx, nSteps);
if (nSteps == MaxSteps)
SetTerminated(idx);
}
else
// If the particle has exited spatial boundary, set corresponding status.
// This is needed for all particle advection cases.
else if (status == ParticleStatus::EXITED_SPATIAL_BOUNDARY)
{
Pos.Set(idx, pt);
SetExitedSpatialBoundary(idx);
}
// If the particle has exited temporal boundary, set corresponding status.
// This is needed when we need temporal interpolation for velocity.
else if (status == ParticleStatus::EXITED_TEMPORAL_BOUNDARY)
{
SetExitedTemporalBoundary(idx);
}
}
/* Set/Change Status */
@ -181,7 +192,8 @@ public:
vtkm::Id GetStatus(const vtkm::Id& idx) const { return Status.Get(idx); }
VTKM_EXEC
T GetTime(const vtkm::Id& idx) const { return Time.Get(idx); }
VTKM_EXEC
void SetTime(const vtkm::Id& idx, T time) const { Time.Set(idx, time); }
protected:
PositionPortal Pos;
@ -232,10 +244,10 @@ public:
Status = statusArray.PrepareForInPlace(DeviceAdapterTag());
Time = timeArray.PrepareForInPlace(DeviceAdapterTag());
MaxSteps = maxSteps;
Length = maxSteps;
vtkm::Id numPos = posArray.GetNumberOfValues();
History = historyArray.PrepareForOutput(numPos * Length, DeviceAdapterTag());
ValidPoint = validPointArray.PrepareForInPlace(DeviceAdapterTag());
Length = maxSteps;
}
VTKM_EXEC_CONT
@ -335,6 +347,8 @@ public:
vtkm::Id GetStatus(const vtkm::Id& idx) const { return Status.Get(idx); }
VTKM_EXEC
T GetTime(const vtkm::Id& idx) const { return Time.Get(idx); }
VTKM_EXEC
void SetTime(const vtkm::Id& idx, T time) const { Time.Set(idx, time); }
vtkm::Vec<T, 3> GetHistory(const vtkm::Id& idx, const vtkm::Id& step) const
{
return History.Get(idx * Length + step);

13
vtkm/worklet/particleadvection/TemporalGridEvaluators.h
View File

@ -56,10 +56,10 @@ public:
throw vtkm::cont::ErrorInternal("Cells are not 3D structured.");
vectors1 = vectorField1.PrepareForInput(DeviceAdapterTag());
vectors2 = vectorField1.PrepareForInput(DeviceAdapterTag());
vectors2 = vectorField2.PrepareForInput(DeviceAdapterTag());
bounds1 = coords1.GetBounds();
bounds2 = coords1.GetBounds();
bounds2 = coords2.GetBounds();
vtkm::cont::CellSetStructured<3> cells1;
vtkm::cont::CellSetStructured<3> cells2;
@ -100,7 +100,7 @@ public:
VTKM_EXEC_CONT
bool IsWithinTemporalBoundary(const FieldType time) const
{
if (time < time1 || time > time2)
if (time < time1 || time >= time2)
return false;
return true;
}
@ -165,9 +165,9 @@ public:
(position[1] - bounds.Y.Min) * scalingFactor[1],
(position[2] - bounds.Z.Min) * scalingFactor[2]);
idx000[0] = floor(scaled[0]);
idx000[1] = floor(scaled[1]);
idx000[2] = floor(scaled[2]);
idx000[0] = floor(scaled[0]) <= dims[0] - 1 ? floor(scaled[0]) : dims[0] - 1;
idx000[1] = floor(scaled[1]) <= dims[1] - 1 ? floor(scaled[1]) : dims[1] - 1;
idx000[2] = floor(scaled[2]) <= dims[2] - 1 ? floor(scaled[2]) : dims[2] - 1;
idx001 = idx000;
idx001[0] = (idx001[0] + 1) <= dims[0] - 1 ? idx001[0] + 1 : dims[0] - 1;
@ -251,6 +251,7 @@ public:
return false;
FieldType proportion = (particleTime - time1) / (time2 - time1);
velocity[0] = (1.0f - proportion) * velocity1[0] + proportion * velocity2[0];
velocity[1] = (1.0f - proportion) * velocity1[1] + proportion * velocity2[1];
velocity[2] = (1.0f - proportion) * velocity1[2] + proportion * velocity2[2];

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

@ -243,7 +243,7 @@ public:
vtkm::worklet::particleadvection::ParticleStatus& status,
vtkm::Vec<FieldType, 3>& pointOut) const
{
status = integrator.Step(pointIn, pointOut);
status = integrator.Step(pointIn, 0, pointOut);
}
private: