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

Cleanup of particle advection worklets

Browse Source 
Created ParticleAdvection.h for worklets that match the worklet pattern.
Slight reorganization of existing code to match this pattern.
Added a unittest particle advection.
This commit is contained in:
Dave Pugmire 2017-07-11 21:22:28 -05:00
parent 77c152abc2
commit 4574f36dec
9 changed files with 328 additions and 77 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
examples/particle_advection/CMakeLists.txt
View File

@ -30,10 +30,10 @@ target_include_directories(Particle_Advection_SERIAL PRIVATE ${VTKm_INCLUDE_DIRS
target_link_libraries(Particle_Advection_SERIAL ${VTKm_LIBRARIES})
target_compile_options(Particle_Advection_SERIAL PRIVATE ${VTKm_COMPILE_OPTIONS})
add_executable(Particle_Advection_RL_SERIAL ParticleAdvectionRectilinear.cxx)
target_include_directories(Particle_Advection_RL_SERIAL PRIVATE ${VTKm_INCLUDE_DIRS})
target_link_libraries(Particle_Advection_RL_SERIAL ${VTKm_LIBRARIES})
target_compile_options(Particle_Advection_RL_SERIAL PRIVATE ${VTKm_COMPILE_OPTIONS})
#add_executable(Particle_Advection_RL_SERIAL ParticleAdvectionRectilinear.cxx)
#target_include_directories(Particle_Advection_RL_SERIAL PRIVATE ${VTKm_INCLUDE_DIRS})
#target_link_libraries(Particle_Advection_RL_SERIAL ${VTKm_LIBRARIES})
#target_compile_options(Particle_Advection_RL_SERIAL PRIVATE ${VTKm_COMPILE_OPTIONS})
if(VTKm_TBB_FOUND)
add_executable(Particle_Advection_TBB ParticleAdvectionTBB.cxx)
@ -41,10 +41,10 @@ if(VTKm_TBB_FOUND)
target_link_libraries(Particle_Advection_TBB ${VTKm_LIBRARIES})
target_compile_options(Particle_Advection_TBB PRIVATE ${VTKm_COMPILE_OPTIONS})
add_executable(Particle_Advection_RL_TBB ParticleAdvectionRectilinearTBB.cxx)
target_include_directories(Particle_Advection_RL_TBB PRIVATE ${VTKm_INCLUDE_DIRS})
target_link_libraries(Particle_Advection_RL_TBB ${VTKm_LIBRARIES})
target_compile_options(Particle_Advection_RL_TBB PRIVATE ${VTKm_COMPILE_OPTIONS})
# add_executable(Particle_Advection_RL_TBB ParticleAdvectionRectilinearTBB.cxx)
# target_include_directories(Particle_Advection_RL_TBB PRIVATE ${VTKm_INCLUDE_DIRS})
# target_link_libraries(Particle_Advection_RL_TBB ${VTKm_LIBRARIES})
# target_compile_options(Particle_Advection_RL_TBB PRIVATE ${VTKm_COMPILE_OPTIONS})
endif()
if(VTKm_CUDA_FOUND)
@ -55,9 +55,9 @@ if(VTKm_CUDA_FOUND)
target_link_libraries(Particle_Advection_CUDA PRIVATE ${VTKm_LIBRARIES})
target_compile_options(Particle_Advection_CUDA PRIVATE ${VTKm_COMPILE_OPTIONS})
cuda_include_directories(${VTKm_INCLUDE_DIRS})
cuda_add_executable(Particle_Advection_RL_CUDA ParticleAdvectionRectilinear.cu)
target_include_directories(Particle_Advection_RL_SERIAL PRIVATE ${VTKm_INCLUDE_DIRS})
target_link_libraries(Particle_Advection_RL_CUDA PRIVATE ${VTKm_LIBRARIES})
target_compile_options(Particle_Advection_RL_CUDA PRIVATE ${VTKm_COMPILE_OPTIONS})
# cuda_include_directories(${VTKm_INCLUDE_DIRS})
# cuda_add_executable(Particle_Advection_RL_CUDA ParticleAdvectionRectilinear.cu)
# target_include_directories(Particle_Advection_RL_SERIAL PRIVATE ${VTKm_INCLUDE_DIRS})
# target_link_libraries(Particle_Advection_RL_CUDA PRIVATE ${VTKm_LIBRARIES})
# target_compile_options(Particle_Advection_RL_CUDA PRIVATE ${VTKm_COMPILE_OPTIONS})
endif()

19
examples/particle_advection/ParticleAdvection.cxx
View File

@ -23,6 +23,7 @@
#endif
#include <vtkm/cont/DataSet.h>
#include <vtkm/worklet/ParticleAdvection.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h>
@ -163,18 +164,24 @@ void RunTest(const std::string& fname,
//time only the actual run
auto t0 = std::chrono::high_resolution_clock::now();
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> seedArray;
seedArray = vtkm::cont::make_ArrayHandle(seeds);
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> fieldArray;
ds.GetField(0).GetData().CopyTo(fieldArray);
if (advectType == 0)
{
vtkm::worklet::particleadvection::ParticleAdvectionWorklet<RK4RGType, FieldType, DeviceAdapter>
pa(rk4, seeds, ds, numSteps);
pa.run(dumpOutput);
vtkm::worklet::ParticleAdvection particleAdvection;
particleAdvection.Run(rk4, seedArray, fieldArray, numSteps, -1, DeviceAdapter());
}
else
{
vtkm::worklet::particleadvection::StreamlineWorklet<RK4RGType, FieldType, DeviceAdapter> sl(
rk4, seeds, ds, numSteps, stepsPerRound, particlesPerRound);
sl.run(dumpOutput);
vtkm::worklet::Streamline streamline;
streamline.Run(
rk4, seedArray, fieldArray, numSteps, stepsPerRound, particlesPerRound, DeviceAdapter());
}
auto t1 = std::chrono::high_resolution_clock::now() - t0;
auto runtime = std::chrono::duration_cast<std::chrono::milliseconds>(t1).count();
std::cerr << "Runtime = " << runtime << " ms " << std::endl;

1
vtkm/worklet/CMakeLists.txt
View File

@ -43,6 +43,7 @@ set(headers
MarchingCubesDataTables.h
Mask.h
MaskPoints.h
ParticleAdvection.h
PointAverage.h
PointElevation.h
RemoveUnusedPoints.h

87
vtkm/worklet/ParticleAdvection.h Normal file
View File

@ -0,0 +1,87 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_ParticleAdvection_h
#define vtk_m_worklet_ParticleAdvection_h
#include <vtkm/Types.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h>
namespace vtkm
{
namespace worklet
{
class ParticleAdvection
{
public:
ParticleAdvection() {}
template <typename IntegratorType,
typename FieldType,
typename PointStorage,
typename FieldStorage,
typename DeviceAdapterTag>
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage> Run(
const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, FieldStorage> fieldArray,
const vtkm::Id& nSteps,
const vtkm::Id& particlesPerRound,
DeviceAdapterTag device)
{
vtkm::worklet::particleadvection::ParticleAdvectionWorklet<IntegratorType,
FieldType,
DeviceAdapterTag>
worklet;
return worklet.Run(it, pts, fieldArray, nSteps, particlesPerRound);
}
};
class Streamline
{
public:
Streamline() {}
template <typename IntegratorType,
typename FieldType,
typename PointStorage,
typename FieldStorage,
typename DeviceAdapterTag>
void Run(const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, FieldStorage> fieldArray,
const vtkm::Id& nSteps,
const vtkm::Id& stepsPerRound,
const vtkm::Id& particlesPerRound,
DeviceAdapterTag device)
{
vtkm::worklet::particleadvection::StreamlineWorklet<IntegratorType, FieldType, DeviceAdapterTag>
worklet;
worklet.Run(it, pts, fieldArray, nSteps, stepsPerRound, particlesPerRound);
}
};
}
}
#endif // vtk_m_worklet_ParticleAdvection_h

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

@ -117,6 +117,9 @@ template <typename PortalType, typename DeviceAdapter, typename FieldType>
class RegularGridEvaluate
{
public:
VTKM_CONT
RegularGridEvaluate() {}
VTKM_CONT
RegularGridEvaluate(const vtkm::cont::DataSet& ds)
{

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

@ -34,6 +34,13 @@ template <typename FieldEvaluateType, typename FieldType, typename PortalType>
class RK4Integrator
{
public:
VTKM_EXEC_CONT
RK4Integrator()
: h(0)
, h_2(0)
{
}
VTKM_EXEC_CONT
RK4Integrator(const FieldEvaluateType& field, FieldType _h)
: f(field)

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

@ -57,7 +57,6 @@ public:
{
vtkm::Vec<FieldType, 3> p = ic.GetPos(idx);
vtkm::Vec<FieldType, 3> p2;
//vtkm::Vec<FieldType, 3> p0 = p;
while (!ic.Done(idx))
{
@ -71,10 +70,6 @@ public:
ic.SetStatusOutOfBounds(idx);
}
}
//std::cerr<<idx<<" DONE"<<std::endl;
//p2 = ic.GetPos(idx);
//std::cerr<<"PIC: "<<idx<<" "<<p0<<" --> "<<p2<<" #steps= "<<ic.GetStep(idx)<<std::endl;
}
ParticleAdvectWorklet(const IntegratorType& it, const FieldPortalConstType& f)
@ -87,6 +82,7 @@ public:
FieldPortalConstType field;
};
template <typename IntegratorType, typename FieldType, typename DeviceAdapterTag>
class ParticleAdvectionWorklet
{
@ -99,42 +95,40 @@ public:
DeviceAdapterTag>
ParticleAdvectWorkletType;
ParticleAdvectionWorklet(const IntegratorType& it,
const std::vector<vtkm::Vec<FieldType, 3>>& pts,
const vtkm::cont::DataSet& _ds,
const vtkm::Id& nSteps,
const vtkm::Id& particlesPerRound = -1)
: integrator(it)
, seeds(pts)
, ds(_ds)
, maxSteps(nSteps)
, ParticlesPerRound(particlesPerRound)
ParticleAdvectionWorklet() {}
template <typename PointStorage, typename FieldStorage>
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage> Run(
const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, FieldStorage> fieldArray,
const vtkm::Id& nSteps,
const vtkm::Id& particlesPerRound = -1)
{
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 3>> fieldArray;
ds.GetField(0).GetData().CopyTo(fieldArray);
integrator = it;
seedArray = pts;
maxSteps = nSteps;
ParticlesPerRound = particlesPerRound;
field = fieldArray.PrepareForInput(DeviceAdapterTag());
return run();
}
~ParticleAdvectionWorklet() {}
void run(bool dumpOutput = false)
private:
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> run(bool dumpOutput = false)
{
typedef typename vtkm::worklet::DispatcherMapField<ParticleAdvectWorkletType>
ParticleWorkletDispatchType;
typedef vtkm::worklet::particleadvection::Particles<FieldType, DeviceAdapterTag> ParticleType;
vtkm::Id totNumSeeds = static_cast<vtkm::Id>(seeds.size());
/*bool NeedParticleRounds = false;*/
vtkm::Id totNumSeeds = static_cast<vtkm::Id>(seedArray.GetNumberOfValues());
vtkm::Id numSeeds = totNumSeeds;
if (ParticlesPerRound == -1 || ParticlesPerRound > totNumSeeds)
numSeeds = totNumSeeds;
else
{
numSeeds = ParticlesPerRound;
/* NeedParticleRounds = true;*/
}
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> posArray =
vtkm::cont::make_ArrayHandle(&seeds[0], numSeeds);
std::vector<vtkm::Id> steps(static_cast<size_t>(numSeeds), 0),
status(static_cast<size_t>(numSeeds), ParticleStatus::OK);
vtkm::cont::ArrayHandle<vtkm::Id> stepArray = vtkm::cont::make_ArrayHandle(&steps[0], numSeeds);
@ -142,7 +136,7 @@ public:
vtkm::cont::make_ArrayHandle(&status[0], numSeeds);
vtkm::cont::ArrayHandleIndex idxArray(numSeeds);
ParticleType particles(posArray, stepArray, statusArray, maxSteps);
ParticleType particles(seedArray, stepArray, statusArray, maxSteps);
ParticleAdvectWorkletType particleWorklet(integrator, field);
ParticleWorkletDispatchType particleWorkletDispatch(particleWorklet);
@ -156,11 +150,12 @@ public:
std::cout << p[0] << " " << p[1] << " " << p[2] << std::endl;
}
}
return seedArray;
}
private:
IntegratorType integrator;
std::vector<vtkm::Vec<FieldType, 3>> seeds;
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> seedArray;
vtkm::cont::DataSet ds;
vtkm::Id maxSteps;
vtkm::Id ParticlesPerRound;
@ -180,26 +175,28 @@ public:
DeviceAdapterTag>
ParticleAdvectWorkletType;
StreamlineWorklet(const IntegratorType& it,
const std::vector<vtkm::Vec<FieldType, 3>>& pts,
const vtkm::cont::DataSet& _ds,
const vtkm::Id& nSteps,
const vtkm::Id stepsPerRound = -1,
const vtkm::Id particlesPerRound = -1)
: integrator(it)
, seeds(pts)
, ds(_ds)
, maxSteps(nSteps)
, StepsPerRound(stepsPerRound)
, ParticlesPerRound(particlesPerRound)
StreamlineWorklet() {}
template <typename PointStorage, typename FieldStorage>
void Run(const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, FieldStorage> fieldArray,
const vtkm::Id& nSteps,
const vtkm::Id stepsPerRound = -1,
const vtkm::Id& particlesPerRound = -1)
{
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 3>> fieldArray;
ds.GetField(0).GetData().CopyTo(fieldArray);
integrator = it;
seedArray = pts;
maxSteps = nSteps;
StepsPerRound = stepsPerRound;
ParticlesPerRound = particlesPerRound;
field = fieldArray.PrepareForInput(DeviceAdapterTag());
run(true);
}
~StreamlineWorklet() {}
private:
void run(bool dumpOutput = false)
{
typedef typename vtkm::worklet::DispatcherMapField<ParticleAdvectWorkletType>
@ -210,17 +207,10 @@ public:
DeviceAdapterTag>
StreamlineRoundType;
vtkm::Id totNumSeeds = static_cast<vtkm::Id>(seeds.size());
vtkm::Id totNumSeeds = static_cast<vtkm::Id>(seedArray.GetNumberOfValues());
bool NeedParticleRounds = false;
/*vtkm::Id numSeeds = totNumSeeds;*/
/*if (ParticlesPerRound == -1 || ParticlesPerRound > totNumSeeds)
numSeeds = totNumSeeds;
else*/
if (!(ParticlesPerRound == -1 || ParticlesPerRound > totNumSeeds))
{
/*numSeeds = ParticlesPerRound;*/
NeedParticleRounds = true;
}
ParticleAdvectWorkletType particleWorklet(integrator, field);
ParticleWorkletDispatchType particleWorkletDispatch(particleWorklet);
@ -239,15 +229,13 @@ public:
if (NeedParticleRounds && num > ParticlesPerRound)
num = ParticlesPerRound;
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> posArray =
vtkm::cont::make_ArrayHandle(&seeds[(size_t)particleOffset], num);
std::vector<vtkm::Id> steps((size_t)num, 0), status((size_t)num, ParticleStatus::OK);
vtkm::cont::ArrayHandle<vtkm::Id> stepArray = vtkm::cont::make_ArrayHandle(&steps[0], num);
vtkm::cont::ArrayHandle<vtkm::Id> statusArray =
vtkm::cont::make_ArrayHandle(&status[0], num);
vtkm::cont::ArrayHandleIndex idxArray(num);
StreamlineType streamlines(posArray, stepArray, statusArray, maxSteps);
StreamlineType streamlines(seedArray, stepArray, statusArray, maxSteps);
particleWorkletDispatch.Invoke(idxArray, streamlines);
if (dumpOutput)
@ -281,8 +269,6 @@ public:
if (NeedParticleRounds && num > ParticlesPerRound)
num = ParticlesPerRound;
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> posArray =
vtkm::cont::make_ArrayHandle(&seeds[(size_t)particleOffset], num);
std::vector<vtkm::Id> steps((size_t)num, 0), status((size_t)num, ParticleStatus::OK);
vtkm::cont::ArrayHandle<vtkm::Id> stepArray = vtkm::cont::make_ArrayHandle(&steps[0], num);
vtkm::cont::ArrayHandle<vtkm::Id> statusArray =
@ -301,7 +287,7 @@ public:
}
StreamlineRoundType streamlines(
posArray, stepArray, statusArray, numSteps, StepsPerRound, stepOffset, maxSteps);
seedArray, stepArray, statusArray, numSteps, StepsPerRound, stepOffset, maxSteps);
particleWorkletDispatch.Invoke(idxArray, streamlines);
auto historyPortal = streamlines.historyArray.GetPortalConstControl();
@ -326,9 +312,8 @@ public:
}
}
private:
IntegratorType integrator;
std::vector<vtkm::Vec<FieldType, 3>> seeds;
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> seedArray;
vtkm::cont::DataSet ds;
vtkm::Id maxSteps;
vtkm::Id StepsPerRound, ParticlesPerRound;

1
vtkm/worklet/testing/CMakeLists.txt
View File

@ -35,6 +35,7 @@ set(unit_tests
UnitTestMarchingCubes.cxx
UnitTestMask.cxx
UnitTestMaskPoints.cxx
UnitTestParticleAdvection.cxx
UnitTestPointElevation.cxx
UnitTestPointGradient.cxx
UnitTestRemoveUnusedPoints.cxx

160
vtkm/worklet/testing/UnitTestParticleAdvection.cxx Normal file
View File

@ -0,0 +1,160 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/DataSetBuilderUniform.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/worklet/ParticleAdvection.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/Particles.h>
vtkm::Float32 vecData[125 * 3] = {
-0.00603248f, -0.0966396f, -0.000732792f, 0.000530014f, -0.0986189f, -0.000806706f,
0.00684929f, -0.100098f, -0.000876566f, 0.0129235f, -0.101102f, -0.000942341f,
0.0187515f, -0.101656f, -0.00100401f, 0.0706091f, -0.083023f, -0.00144278f,
0.0736404f, -0.0801616f, -0.00145784f, 0.0765194f, -0.0772063f, -0.00147036f,
0.0792559f, -0.0741751f, -0.00148051f, 0.0818589f, -0.071084f, -0.00148843f,
0.103585f, -0.0342287f, -0.001425f, 0.104472f, -0.0316147f, -0.00140433f,
0.105175f, -0.0291574f, -0.00138057f, 0.105682f, -0.0268808f, -0.00135357f,
0.105985f, -0.0248099f, -0.00132315f, -0.00244603f, -0.0989576f, -0.000821705f,
0.00389525f, -0.100695f, -0.000894513f, 0.00999301f, -0.10193f, -0.000963114f,
0.0158452f, -0.102688f, -0.00102747f, 0.0214509f, -0.102995f, -0.00108757f,
0.0708166f, -0.081799f, -0.00149941f, 0.0736939f, -0.0787879f, -0.00151236f,
0.0764359f, -0.0756944f, -0.00152297f, 0.0790546f, -0.0725352f, -0.00153146f,
0.0815609f, -0.0693255f, -0.001538f, -0.00914287f, -0.104658f, -0.001574f,
-0.00642891f, -0.10239f, -0.00159659f, -0.00402289f, -0.0994835f, -0.00160731f,
-0.00194792f, -0.0959752f, -0.00160528f, -0.00022818f, -0.0919077f, -0.00158957f,
-0.0134913f, -0.0274735f, -9.50056e-05f, -0.0188683f, -0.023273f, 0.000194107f,
-0.0254516f, -0.0197589f, 0.000529693f, -0.0312798f, -0.0179514f, 0.00083619f,
-0.0360426f, -0.0177537f, 0.00110164f, 0.0259929f, -0.0204479f, -0.000304646f,
0.033336f, -0.0157385f, -0.000505569f, 0.0403427f, -0.0104637f, -0.000693529f,
0.0469371f, -0.00477766f, -0.000865609f, 0.0530722f, 0.0011701f, -0.00102f,
-0.0121869f, -0.10317f, -0.0015868f, -0.0096549f, -0.100606f, -0.00160377f,
-0.00743038f, -0.0973796f, -0.00160783f, -0.00553901f, -0.0935261f, -0.00159792f,
-0.00400821f, -0.0890871f, -0.00157287f, -0.0267803f, -0.0165823f, 0.000454173f,
-0.0348303f, -0.011642f, 0.000881271f, -0.0424964f, -0.00870761f, 0.00129226f,
-0.049437f, -0.00781358f, 0.0016728f, -0.0552635f, -0.00888708f, 0.00200659f,
-0.0629746f, -0.0721524f, -0.00160475f, -0.0606813f, -0.0677576f, -0.00158427f,
-0.0582203f, -0.0625009f, -0.00154304f, -0.0555686f, -0.0563905f, -0.00147822f,
-0.0526988f, -0.0494369f, -0.00138643f, 0.0385695f, 0.115704f, 0.00674413f,
0.056434f, 0.128273f, 0.00869052f, 0.0775564f, 0.137275f, 0.0110399f,
0.102515f, 0.140823f, 0.0138637f, 0.131458f, 0.136024f, 0.0171804f,
0.0595175f, -0.0845927f, 0.00512454f, 0.0506615f, -0.0680369f, 0.00376604f,
0.0434904f, -0.0503557f, 0.00261592f, 0.0376711f, -0.0318716f, 0.00163301f,
0.0329454f, -0.0128019f, 0.000785352f, -0.0664062f, -0.0701094f, -0.00160644f,
-0.0641074f, -0.0658893f, -0.00158969f, -0.0616054f, -0.0608302f, -0.00155303f,
-0.0588734f, -0.0549447f, -0.00149385f, -0.0558797f, -0.0482482f, -0.00140906f,
0.0434062f, 0.102969f, 0.00581269f, 0.0619547f, 0.112838f, 0.00742057f,
0.0830229f, 0.118752f, 0.00927516f, 0.106603f, 0.119129f, 0.0113757f,
0.132073f, 0.111946f, 0.0136613f, -0.0135758f, -0.0934604f, -0.000533868f,
-0.00690763f, -0.0958773f, -0.000598878f, -0.000475275f, -0.0977838f, -0.000660985f,
0.00571866f, -0.0992032f, -0.0007201f, 0.0116724f, -0.10016f, -0.000776144f,
0.0651428f, -0.0850475f, -0.00120243f, 0.0682895f, -0.0823666f, -0.00121889f,
0.0712792f, -0.0795772f, -0.00123291f, 0.0741224f, -0.0766981f, -0.00124462f,
0.076829f, -0.0737465f, -0.00125416f, 0.10019f, -0.0375515f, -0.00121866f,
0.101296f, -0.0348723f, -0.00120216f, 0.102235f, -0.0323223f, -0.00118309f,
0.102994f, -0.0299234f, -0.00116131f, 0.103563f, -0.0276989f, -0.0011367f,
-0.00989236f, -0.0958821f, -0.000608883f, -0.00344154f, -0.0980645f, -0.000673641f,
0.00277318f, -0.0997337f, -0.000735354f, 0.00874908f, -0.100914f, -0.000793927f,
0.0144843f, -0.101629f, -0.000849279f, 0.0654428f, -0.0839355f, -0.00125739f,
0.0684225f, -0.0810989f, -0.00127208f, 0.0712599f, -0.0781657f, -0.00128444f,
0.0739678f, -0.0751541f, -0.00129465f, 0.076558f, -0.0720804f, -0.00130286f,
-0.0132841f, -0.103948f, -0.00131159f, -0.010344f, -0.102328f, -0.0013452f,
-0.00768637f, -0.100054f, -0.00136938f, -0.00533293f, -0.0971572f, -0.00138324f,
-0.00330643f, -0.0936735f, -0.00138586f, -0.0116984f, -0.0303752f, -0.000229102f,
-0.0149879f, -0.0265231f, -3.43823e-05f, -0.0212917f, -0.0219544f, 0.000270283f,
-0.0277756f, -0.0186879f, 0.000582781f, -0.0335115f, -0.0171098f, 0.00086919f,
0.0170095f, -0.025299f, -3.73557e-05f, 0.024552f, -0.0214351f, -0.000231975f,
0.0318714f, -0.0168568f, -0.000417463f, 0.0388586f, -0.0117131f, -0.000589883f,
0.0454388f, -0.00615626f, -0.000746594f, -0.0160785f, -0.102675f, -0.00132891f,
-0.0133174f, -0.100785f, -0.00135859f, -0.0108365f, -0.0982184f, -0.00137801f,
-0.00865931f, -0.0950053f, -0.00138614f, -0.00681126f, -0.0911806f, -0.00138185f,
-0.0208973f, -0.0216631f, 0.000111231f, -0.0289373f, -0.0151081f, 0.000512553f,
-0.0368736f, -0.0104306f, 0.000911793f, -0.0444294f, -0.00773838f, 0.00129762f,
-0.0512663f, -0.00706554f, 0.00165611f
};
void TestParticleAdvection()
{
typedef VTKM_DEFAULT_DEVICE_ADAPTER_TAG DeviceAdapter;
typedef vtkm::Float32 FieldType;
typedef vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> FieldHandle;
typedef
typename FieldHandle::template ExecutionTypes<DeviceAdapter>::PortalConst FieldPortalConstType;
std::cout << "Testing Integrators for ParticleAdvection Worklet" << std::endl;
FieldType stepSize = 0.01f;
vtkm::cont::DataSetBuilderUniform dataSetBuilder;
const vtkm::Id3 dims(5, 5, 5);
vtkm::Id nElements = dims[0] * dims[1] * dims[2] * 3;
std::vector<vtkm::Vec<vtkm::Float32, 3>> field;
for (vtkm::Id i = 0; i < nElements; i++)
{
FieldType x = vecData[i];
FieldType y = vecData[++i];
FieldType z = vecData[++i];
vtkm::Vec<FieldType, 3> vec(x, y, z);
field.push_back(vtkm::Normal(vec));
}
vtkm::cont::DataSet ds = dataSetBuilder.Create(dims);
vtkm::worklet::particleadvection::RegularGridEvaluate<FieldPortalConstType,
DeviceAdapter,
FieldType>
eval(ds);
typedef vtkm::worklet::particleadvection::RegularGridEvaluate<FieldPortalConstType,
DeviceAdapter,
FieldType>
RGEvalType;
typedef vtkm::worklet::particleadvection::RK4Integrator<RGEvalType,
FieldType,
FieldPortalConstType>
RK4RGType;
RK4RGType rk4(eval, stepSize);
std::vector<vtkm::Vec<FieldType, 3>> pts;
pts.push_back(vtkm::Vec<FieldType, 3>(1, 1, 1));
pts.push_back(vtkm::Vec<FieldType, 3>(2, 2, 2));
pts.push_back(vtkm::Vec<FieldType, 3>(3, 3, 3));
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> seeds;
seeds = vtkm::cont::make_ArrayHandle(pts);
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> fieldArray;
fieldArray = vtkm::cont::make_ArrayHandle(field);
vtkm::worklet::ParticleAdvection particleAdvection;
auto endSeeds = particleAdvection.Run(rk4, seeds, fieldArray, 100, -1, DeviceAdapter());
VTKM_TEST_ASSERT(seeds.GetNumberOfValues() == endSeeds.GetNumberOfValues(),
"Number of output particles does not match input.");
}
int UnitTestParticleAdvection(int, char* [])
{
return vtkm::cont::testing::Testing::Run(TestParticleAdvection);
}