Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-08 13:23:51 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Fix for small euler step for particle advection.

Browse Source
This commit is contained in:
dpugmire 2020-04-23 16:14:28 -04:00
parent d22b658c3a
commit d2e9b3d306
5 changed files with 84 additions and 75 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
vtkm/exec/CellLocatorRectilinearGrid.h
View File

@ -117,6 +117,7 @@ public:
if (point[dim] == MaxPoint[dim])
{
logicalCell[dim] = this->PointDimensions[dim] - 2;
parametric[dim] = static_cast<vtkm::FloatDefault>(1);
continue;
}

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

@ -112,6 +112,7 @@ public:
for (vtkm::IdComponent i = 0; i < nVerts; i++)
fieldValues.Append(Field.Get(ptIndices[i]));
vtkm::exec::CellInterpolate(fieldValues, parametric, cellShape, out);
status.SetOk();

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

@ -13,7 +13,6 @@
#ifndef vtk_m_worklet_particleadvection_Integrators_h
#define vtk_m_worklet_particleadvection_Integrators_h
#include <iomanip>
#include <limits>
#include <vtkm/Bitset.h>
@ -168,43 +167,57 @@ protected:
//we terminate, the outside value is the stepsize that will nudge
//the particle outside the dataset.
//The binary search will be between [0, this->StepLength]
vtkm::FloatDefault stepShort = 0, stepLong = this->StepLength;
vtkm::Vec3f currPos(inpos), velocity;
vtkm::FloatDefault currTime = time;
//The binary search will be between {0, this->StepLength}
vtkm::FloatDefault stepRange[2] = { 0, this->StepLength };
auto evalStatus = this->Evaluator.Evaluate(currPos, time, velocity);
vtkm::Vec3f currPos(inpos), currVel;
auto evalStatus = this->Evaluator.Evaluate(currPos, time, currVel);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
const vtkm::FloatDefault eps = vtkm::Epsilon<vtkm::FloatDefault>();
const vtkm::FloatDefault eps = vtkm::Epsilon<vtkm::FloatDefault>() * 10;
vtkm::FloatDefault div = 1;
for (int i = 0; i < 50; i++)
vtkm::Vec3f tmp;
while ((stepRange[1] - stepRange[0]) > eps)
{
//Try a step midway between stepRange[0] and stepRange[1]
div *= 2;
vtkm::FloatDefault stepCurr = stepShort + (this->StepLength / div);
vtkm::FloatDefault stepCurr = stepRange[0] + (this->StepLength / div);
//See if we can step by stepCurr.
IntegratorStatus status = this->CheckStep(inpos, stepCurr, time, velocity);
if (status.CheckOk())
IntegratorStatus status = this->CheckStep(inpos, stepCurr, time, currVel);
if (status.CheckOk()) //Integration step succedded.
{
currPos = inpos + velocity * stepShort;
stepShort = stepCurr;
//See if this point is in/out.
auto newPos = inpos + stepCurr * currVel;
evalStatus = this->Evaluator.Evaluate(newPos, time + stepCurr, tmp);
if (evalStatus.CheckOk())
{
//Point still in. Update currPos and set range to {stepCurr, stepRange[1]}
currPos = newPos;
stepRange[0] = stepCurr;
}
else
{
//The step succedded, but the next point is outside.
//Step too long. Set range to: {stepRange[0], stepCurr} and continue.
stepRange[1] = stepCurr;
}
}
else
stepLong = stepCurr;
//Stop if step bracket is small enough.
if (stepLong - stepShort < eps)
break;
{
//Step too long. Set range to: {stepRange[0], stepCurr} and continue.
stepRange[1] = stepCurr;
}
}
//Take Euler step.
currTime = time + stepShort;
evalStatus = this->Evaluator.Evaluate(currPos, currTime, velocity);
evalStatus = this->Evaluator.Evaluate(currPos, time + stepRange[0], currVel);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
outpos = currPos + stepLong * velocity;
//Update the position and time.
outpos = currPos + stepRange[1] * currVel;
time += stepRange[1];
return IntegratorStatus(true, true, !this->Evaluator.IsWithinTemporalBoundary(time));
}
@ -292,6 +305,12 @@ public:
if ((time + stepLength + vtkm::Epsilon<vtkm::FloatDefault>() - boundary) > 0.0)
stepLength = boundary - time;
//k1 = F(p,t)
//k2 = F(p+hk1/2, t+h/2
//k3 = F(p+hk2/2, t+h/2
//k4 = F(p+hk3, t+h)
//Yn+1 = Yn + 1/6 h (k1+2k2+2k3+k4)
vtkm::FloatDefault var1 = (stepLength / static_cast<vtkm::FloatDefault>(2));
vtkm::FloatDefault var2 = time + var1;
vtkm::FloatDefault var3 = time + stepLength;
@ -313,7 +332,7 @@ public:
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
velocity = (k1 + 2 * k2 + 2 * k3 + k4) / 6.0f;
velocity = (k1 + 2 * k2 + 2 * k3 + k4) / static_cast<vtkm::FloatDefault>(6);
return IntegratorStatus(true, false, evalStatus.CheckTemporalBounds());
}
};

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

@ -51,7 +51,7 @@ public:
{
vtkm::Particle particle = integralCurve.GetParticle(idx);
vtkm::Vec3f inpos = particle.Pos;
vtkm::Vec3f currPos = particle.Pos;
vtkm::FloatDefault time = particle.Time;
bool tookAnySteps = false;
@ -63,26 +63,25 @@ public:
integralCurve.PreStepUpdate(idx);
do
{
//vtkm::Particle p = integralCurve.GetParticle(idx);
//std::cout<<idx<<": "<<inpos<<" #"<<p.NumSteps<<" "<<p.Status<<std::endl;
vtkm::Vec3f outpos;
auto status = integrator->Step(inpos, time, outpos);
auto status = integrator->Step(currPos, time, outpos);
if (status.CheckOk())
{
integralCurve.StepUpdate(idx, time, outpos);
tookAnySteps = true;
inpos = outpos;
currPos = outpos;
}
//We can't take a step inside spatial boundary.
//Try and take a step just past the boundary.
else if (status.CheckSpatialBounds())
{
IntegratorStatus status2 = integrator->SmallStep(inpos, time, outpos);
IntegratorStatus status2 = integrator->SmallStep(currPos, time, outpos);
if (status2.CheckOk())
{
integralCurve.StepUpdate(idx, time, outpos);
tookAnySteps = true;
currPos = outpos;
//we took a step, so use this status to consider below.
status = status2;
@ -95,9 +94,6 @@ public:
//Mark if any steps taken
integralCurve.UpdateTookSteps(idx, tookAnySteps);
//particle = integralCurve.GetParticle(idx);
//std::cout<<idx<<": "<<inpos<<" #"<<particle.NumSteps<<" "<<particle.Status<<std::endl;
}
};

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

@ -21,13 +21,7 @@
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/Particles.h>
#include <vtkm/io/writer/VTKDataSetWriter.h>
//timers
#include <chrono>
#include <ctime>
#include <random>
namespace
{
@ -98,20 +92,29 @@ vtkm::FloatDefault vecData[125 * 3] = {
};
}
vtkm::Vec3f RandomPoint(const vtkm::Bounds& bounds)
void GenerateRandomParticles(std::vector<vtkm::Particle>& points,
std::size_t N,
vtkm::Bounds& bounds,
std::size_t seed = 314)
{
vtkm::FloatDefault rx =
static_cast<vtkm::FloatDefault>(rand()) / static_cast<vtkm::FloatDefault>(RAND_MAX);
vtkm::FloatDefault ry =
static_cast<vtkm::FloatDefault>(rand()) / static_cast<vtkm::FloatDefault>(RAND_MAX);
vtkm::FloatDefault rz =
static_cast<vtkm::FloatDefault>(rand()) / static_cast<vtkm::FloatDefault>(RAND_MAX);
std::random_device device;
std::default_random_engine generator(static_cast<vtkm::UInt32>(seed));
vtkm::FloatDefault zero(0), one(1);
std::uniform_real_distribution<vtkm::FloatDefault> distribution(zero, one);
vtkm::Vec3f p;
p[0] = static_cast<vtkm::FloatDefault>(bounds.X.Min + rx * bounds.X.Length());
p[1] = static_cast<vtkm::FloatDefault>(bounds.Y.Min + ry * bounds.Y.Length());
p[2] = static_cast<vtkm::FloatDefault>(bounds.Z.Min + rz * bounds.Z.Length());
return p;
points.resize(0);
for (std::size_t i = 0; i < N; i++)
{
vtkm::FloatDefault rx = distribution(generator);
vtkm::FloatDefault ry = distribution(generator);
vtkm::FloatDefault rz = distribution(generator);
vtkm::Vec3f p;
p[0] = static_cast<vtkm::FloatDefault>(bounds.X.Min + rx * bounds.X.Length());
p[1] = static_cast<vtkm::FloatDefault>(bounds.Y.Min + ry * bounds.Y.Length());
p[2] = static_cast<vtkm::FloatDefault>(bounds.Z.Min + rz * bounds.Z.Length());
points.push_back(vtkm::Particle(p, static_cast<vtkm::Id>(i)));
}
}
vtkm::cont::DataSet CreateUniformDataSet(const vtkm::Bounds& bounds, const vtkm::Id3& dims)
@ -375,7 +378,7 @@ void ValidateIntegrator(const IntegratorType& integrator,
VTKM_TEST_ASSERT(result == pointsPortal.Get(index).Pos,
"Error in evaluator result for [OUTSIDE SPATIAL]" + msg);
else
VTKM_TEST_ASSERT(result == expStepResults[(size_t)index],
VTKM_TEST_ASSERT(result == expStepResults[static_cast<size_t>(index)],
"Error in evaluator result for " + msg);
}
}
@ -401,16 +404,13 @@ void ValidateIntegratorForBoundary(const vtkm::Bounds& bounds,
for (vtkm::Id index = 0; index < numPoints; index++)
{
Status status = statusPortal.Get(index);
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.
vtkm::Vec3f result = resultsPortal.Get(index);
if (bounds.Contains(result))
{
std::cout << "Failure. " << bounds << std::endl;
std::cout << std::setprecision(12) << index << ": " << bounds << " res= " << result
<< std::endl;
}
//VTKM_TEST_ASSERT(!bounds.Contains(result), "Tolerance not satisfied.");
VTKM_TEST_ASSERT(!bounds.Contains(result),
"Integrator did not step out of boundary for " + msg);
}
}
@ -438,7 +438,6 @@ void TestEvaluators()
std::vector<vtkm::Id3> dims;
dims.push_back(vtkm::Id3(5, 5, 5));
dims.push_back(vtkm::Id3(10, 5, 5));
dims.push_back(vtkm::Id3(10, 5, 5));
for (auto& dim : dims)
{
@ -457,8 +456,7 @@ void TestEvaluators()
vtkm::FloatDefault stepSize = 0.1f;
std::vector<vtkm::Particle> pointIns;
std::vector<vtkm::Vec3f> stepResult;
//Create a bunch of random points in the bounds.
srand(314);
//Generate points 2 steps inside the bounding box.
vtkm::Bounds interiorBounds = bound;
interiorBounds.X.Min += 2 * stepSize;
@ -467,12 +465,10 @@ void TestEvaluators()
interiorBounds.X.Max -= 2 * stepSize;
interiorBounds.Y.Max -= 2 * stepSize;
interiorBounds.Z.Max -= 2 * stepSize;
for (int k = 0; k < 38; k++)
{
auto p = RandomPoint(interiorBounds);
pointIns.push_back(vtkm::Particle(p, k));
stepResult.push_back(p + vec * stepSize);
}
GenerateRandomParticles(pointIns, 38, interiorBounds);
for (auto& p : pointIns)
stepResult.push_back(p.Pos + vec * stepSize);
vtkm::Range xRange, yRange, zRange;
@ -496,8 +492,7 @@ void TestEvaluators()
// All velocities are in the +ve direction.
std::vector<vtkm::Particle> boundaryPoints;
for (int k = 0; k < 10; k++)
boundaryPoints.push_back(vtkm::Particle(RandomPoint(forBoundary), k));
GenerateRandomParticles(boundaryPoints, 10, forBoundary, 919);
for (auto& ds : dataSets)
{
@ -506,7 +501,6 @@ void TestEvaluators()
RK4Type rk4(gridEval, stepSize);
ValidateIntegrator(rk4, pointIns, stepResult, "constant vector RK4");
ValidateIntegratorForBoundary(bound, rk4, boundaryPoints, "constant vector RK4");
}
}
@ -588,9 +582,7 @@ void TestParticleWorkletsWithDataSetTypes()
//Generate three random points.
std::vector<vtkm::Particle> pts;
pts.push_back(vtkm::Particle(RandomPoint(bound), 0));
pts.push_back(vtkm::Particle(RandomPoint(bound), 1));
pts.push_back(vtkm::Particle(RandomPoint(bound), 2));
GenerateRandomParticles(pts, 3, bound, 111);
std::vector<vtkm::Particle> pts2 = pts;
vtkm::Id nSeeds = static_cast<vtkm::Id>(pts.size());