mirror of
https://gitlab.kitware.com/vtk/vtk-m
synced 2024-09-16 17:22:55 +00:00
Merge topic 'fix_euler_step_particleAdvection2'
2ecca9edf Merge branch 'master' of https://gitlab.kitware.com/vtk/vtk-m into fix_euler_step_particleAdvection2 d2e9b3d30 Fix for small euler step for particle advection. Acked-by: Kitware Robot <kwrobot@kitware.com> Merge-request: !2057
This commit is contained in:
commit
87470a89b7
@ -117,6 +117,7 @@ public:
|
|||||||
if (point[dim] == MaxPoint[dim])
|
if (point[dim] == MaxPoint[dim])
|
||||||
{
|
{
|
||||||
logicalCell[dim] = this->PointDimensions[dim] - 2;
|
logicalCell[dim] = this->PointDimensions[dim] - 2;
|
||||||
|
parametric[dim] = static_cast<vtkm::FloatDefault>(1);
|
||||||
continue;
|
continue;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -112,6 +112,7 @@ public:
|
|||||||
|
|
||||||
for (vtkm::IdComponent i = 0; i < nVerts; i++)
|
for (vtkm::IdComponent i = 0; i < nVerts; i++)
|
||||||
fieldValues.Append(Field.Get(ptIndices[i]));
|
fieldValues.Append(Field.Get(ptIndices[i]));
|
||||||
|
|
||||||
vtkm::exec::CellInterpolate(fieldValues, parametric, cellShape, out);
|
vtkm::exec::CellInterpolate(fieldValues, parametric, cellShape, out);
|
||||||
|
|
||||||
status.SetOk();
|
status.SetOk();
|
||||||
|
@ -13,7 +13,6 @@
|
|||||||
#ifndef vtk_m_worklet_particleadvection_Integrators_h
|
#ifndef vtk_m_worklet_particleadvection_Integrators_h
|
||||||
#define vtk_m_worklet_particleadvection_Integrators_h
|
#define vtk_m_worklet_particleadvection_Integrators_h
|
||||||
|
|
||||||
#include <iomanip>
|
|
||||||
#include <limits>
|
#include <limits>
|
||||||
|
|
||||||
#include <vtkm/Bitset.h>
|
#include <vtkm/Bitset.h>
|
||||||
@ -168,43 +167,57 @@ protected:
|
|||||||
//we terminate, the outside value is the stepsize that will nudge
|
//we terminate, the outside value is the stepsize that will nudge
|
||||||
//the particle outside the dataset.
|
//the particle outside the dataset.
|
||||||
|
|
||||||
//The binary search will be between [0, this->StepLength]
|
//The binary search will be between {0, this->StepLength}
|
||||||
vtkm::FloatDefault stepShort = 0, stepLong = this->StepLength;
|
vtkm::FloatDefault stepRange[2] = { 0, this->StepLength };
|
||||||
vtkm::Vec3f currPos(inpos), velocity;
|
|
||||||
vtkm::FloatDefault currTime = time;
|
|
||||||
|
|
||||||
auto evalStatus = this->Evaluator.Evaluate(currPos, time, velocity);
|
vtkm::Vec3f currPos(inpos), currVel;
|
||||||
|
auto evalStatus = this->Evaluator.Evaluate(currPos, time, currVel);
|
||||||
if (evalStatus.CheckFail())
|
if (evalStatus.CheckFail())
|
||||||
return IntegratorStatus(evalStatus);
|
return IntegratorStatus(evalStatus);
|
||||||
|
|
||||||
const vtkm::FloatDefault eps = vtkm::Epsilon<vtkm::FloatDefault>();
|
const vtkm::FloatDefault eps = vtkm::Epsilon<vtkm::FloatDefault>() * 10;
|
||||||
vtkm::FloatDefault div = 1;
|
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;
|
div *= 2;
|
||||||
vtkm::FloatDefault stepCurr = stepShort + (this->StepLength / div);
|
vtkm::FloatDefault stepCurr = stepRange[0] + (this->StepLength / div);
|
||||||
|
|
||||||
//See if we can step by stepCurr.
|
//See if we can step by stepCurr.
|
||||||
IntegratorStatus status = this->CheckStep(inpos, stepCurr, time, velocity);
|
IntegratorStatus status = this->CheckStep(inpos, stepCurr, time, currVel);
|
||||||
if (status.CheckOk())
|
|
||||||
|
if (status.CheckOk()) //Integration step succedded.
|
||||||
{
|
{
|
||||||
currPos = inpos + velocity * stepShort;
|
//See if this point is in/out.
|
||||||
stepShort = stepCurr;
|
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
|
else
|
||||||
stepLong = stepCurr;
|
{
|
||||||
|
//Step too long. Set range to: {stepRange[0], stepCurr} and continue.
|
||||||
//Stop if step bracket is small enough.
|
stepRange[1] = stepCurr;
|
||||||
if (stepLong - stepShort < eps)
|
}
|
||||||
break;
|
|
||||||
}
|
}
|
||||||
|
evalStatus = this->Evaluator.Evaluate(currPos, time + stepRange[0], currVel);
|
||||||
//Take Euler step.
|
|
||||||
currTime = time + stepShort;
|
|
||||||
evalStatus = this->Evaluator.Evaluate(currPos, currTime, velocity);
|
|
||||||
if (evalStatus.CheckFail())
|
if (evalStatus.CheckFail())
|
||||||
return IntegratorStatus(evalStatus);
|
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));
|
return IntegratorStatus(true, true, !this->Evaluator.IsWithinTemporalBoundary(time));
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -292,6 +305,12 @@ public:
|
|||||||
if ((time + stepLength + vtkm::Epsilon<vtkm::FloatDefault>() - boundary) > 0.0)
|
if ((time + stepLength + vtkm::Epsilon<vtkm::FloatDefault>() - boundary) > 0.0)
|
||||||
stepLength = boundary - time;
|
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 var1 = (stepLength / static_cast<vtkm::FloatDefault>(2));
|
||||||
vtkm::FloatDefault var2 = time + var1;
|
vtkm::FloatDefault var2 = time + var1;
|
||||||
vtkm::FloatDefault var3 = time + stepLength;
|
vtkm::FloatDefault var3 = time + stepLength;
|
||||||
@ -313,7 +332,7 @@ public:
|
|||||||
if (evalStatus.CheckFail())
|
if (evalStatus.CheckFail())
|
||||||
return IntegratorStatus(evalStatus);
|
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());
|
return IntegratorStatus(true, false, evalStatus.CheckTemporalBounds());
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
|
@ -51,7 +51,7 @@ public:
|
|||||||
{
|
{
|
||||||
vtkm::Particle particle = integralCurve.GetParticle(idx);
|
vtkm::Particle particle = integralCurve.GetParticle(idx);
|
||||||
|
|
||||||
vtkm::Vec3f inpos = particle.Pos;
|
vtkm::Vec3f currPos = particle.Pos;
|
||||||
vtkm::FloatDefault time = particle.Time;
|
vtkm::FloatDefault time = particle.Time;
|
||||||
bool tookAnySteps = false;
|
bool tookAnySteps = false;
|
||||||
|
|
||||||
@ -63,26 +63,25 @@ public:
|
|||||||
integralCurve.PreStepUpdate(idx);
|
integralCurve.PreStepUpdate(idx);
|
||||||
do
|
do
|
||||||
{
|
{
|
||||||
//vtkm::Particle p = integralCurve.GetParticle(idx);
|
|
||||||
//std::cout<<idx<<": "<<inpos<<" #"<<p.NumSteps<<" "<<p.Status<<std::endl;
|
|
||||||
vtkm::Vec3f outpos;
|
vtkm::Vec3f outpos;
|
||||||
auto status = integrator->Step(inpos, time, outpos);
|
auto status = integrator->Step(currPos, time, outpos);
|
||||||
if (status.CheckOk())
|
if (status.CheckOk())
|
||||||
{
|
{
|
||||||
integralCurve.StepUpdate(idx, time, outpos);
|
integralCurve.StepUpdate(idx, time, outpos);
|
||||||
tookAnySteps = true;
|
tookAnySteps = true;
|
||||||
inpos = outpos;
|
currPos = outpos;
|
||||||
}
|
}
|
||||||
|
|
||||||
//We can't take a step inside spatial boundary.
|
//We can't take a step inside spatial boundary.
|
||||||
//Try and take a step just past the boundary.
|
//Try and take a step just past the boundary.
|
||||||
else if (status.CheckSpatialBounds())
|
else if (status.CheckSpatialBounds())
|
||||||
{
|
{
|
||||||
IntegratorStatus status2 = integrator->SmallStep(inpos, time, outpos);
|
IntegratorStatus status2 = integrator->SmallStep(currPos, time, outpos);
|
||||||
if (status2.CheckOk())
|
if (status2.CheckOk())
|
||||||
{
|
{
|
||||||
integralCurve.StepUpdate(idx, time, outpos);
|
integralCurve.StepUpdate(idx, time, outpos);
|
||||||
tookAnySteps = true;
|
tookAnySteps = true;
|
||||||
|
currPos = outpos;
|
||||||
|
|
||||||
//we took a step, so use this status to consider below.
|
//we took a step, so use this status to consider below.
|
||||||
status = status2;
|
status = status2;
|
||||||
@ -95,9 +94,6 @@ public:
|
|||||||
|
|
||||||
//Mark if any steps taken
|
//Mark if any steps taken
|
||||||
integralCurve.UpdateTookSteps(idx, tookAnySteps);
|
integralCurve.UpdateTookSteps(idx, tookAnySteps);
|
||||||
|
|
||||||
//particle = integralCurve.GetParticle(idx);
|
|
||||||
//std::cout<<idx<<": "<<inpos<<" #"<<particle.NumSteps<<" "<<particle.Status<<std::endl;
|
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
|
|
||||||
|
@ -21,13 +21,7 @@
|
|||||||
#include <vtkm/worklet/particleadvection/Integrators.h>
|
#include <vtkm/worklet/particleadvection/Integrators.h>
|
||||||
#include <vtkm/worklet/particleadvection/Particles.h>
|
#include <vtkm/worklet/particleadvection/Particles.h>
|
||||||
|
|
||||||
#include <vtkm/io/writer/VTKDataSetWriter.h>
|
#include <random>
|
||||||
|
|
||||||
//timers
|
|
||||||
#include <chrono>
|
|
||||||
#include <ctime>
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
namespace
|
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 =
|
std::random_device device;
|
||||||
static_cast<vtkm::FloatDefault>(rand()) / static_cast<vtkm::FloatDefault>(RAND_MAX);
|
std::default_random_engine generator(static_cast<vtkm::UInt32>(seed));
|
||||||
vtkm::FloatDefault ry =
|
vtkm::FloatDefault zero(0), one(1);
|
||||||
static_cast<vtkm::FloatDefault>(rand()) / static_cast<vtkm::FloatDefault>(RAND_MAX);
|
std::uniform_real_distribution<vtkm::FloatDefault> distribution(zero, one);
|
||||||
vtkm::FloatDefault rz =
|
|
||||||
static_cast<vtkm::FloatDefault>(rand()) / static_cast<vtkm::FloatDefault>(RAND_MAX);
|
|
||||||
|
|
||||||
vtkm::Vec3f p;
|
points.resize(0);
|
||||||
p[0] = static_cast<vtkm::FloatDefault>(bounds.X.Min + rx * bounds.X.Length());
|
for (std::size_t i = 0; i < N; i++)
|
||||||
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());
|
vtkm::FloatDefault rx = distribution(generator);
|
||||||
return p;
|
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)
|
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,
|
VTKM_TEST_ASSERT(result == pointsPortal.Get(index).Pos,
|
||||||
"Error in evaluator result for [OUTSIDE SPATIAL]" + msg);
|
"Error in evaluator result for [OUTSIDE SPATIAL]" + msg);
|
||||||
else
|
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);
|
"Error in evaluator result for " + msg);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
@ -401,16 +404,13 @@ void ValidateIntegratorForBoundary(const vtkm::Bounds& bounds,
|
|||||||
for (vtkm::Id index = 0; index < numPoints; index++)
|
for (vtkm::Id index = 0; index < numPoints; index++)
|
||||||
{
|
{
|
||||||
Status status = statusPortal.Get(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);
|
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);
|
vtkm::Vec3f result = resultsPortal.Get(index);
|
||||||
if (bounds.Contains(result))
|
VTKM_TEST_ASSERT(!bounds.Contains(result),
|
||||||
{
|
"Integrator did not step out of boundary for " + msg);
|
||||||
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.");
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -438,7 +438,6 @@ void TestEvaluators()
|
|||||||
std::vector<vtkm::Id3> dims;
|
std::vector<vtkm::Id3> dims;
|
||||||
dims.push_back(vtkm::Id3(5, 5, 5));
|
dims.push_back(vtkm::Id3(5, 5, 5));
|
||||||
dims.push_back(vtkm::Id3(10, 5, 5));
|
dims.push_back(vtkm::Id3(10, 5, 5));
|
||||||
dims.push_back(vtkm::Id3(10, 5, 5));
|
|
||||||
|
|
||||||
for (auto& dim : dims)
|
for (auto& dim : dims)
|
||||||
{
|
{
|
||||||
@ -457,8 +456,7 @@ void TestEvaluators()
|
|||||||
vtkm::FloatDefault stepSize = 0.1f;
|
vtkm::FloatDefault stepSize = 0.1f;
|
||||||
std::vector<vtkm::Particle> pointIns;
|
std::vector<vtkm::Particle> pointIns;
|
||||||
std::vector<vtkm::Vec3f> stepResult;
|
std::vector<vtkm::Vec3f> stepResult;
|
||||||
//Create a bunch of random points in the bounds.
|
|
||||||
srand(314);
|
|
||||||
//Generate points 2 steps inside the bounding box.
|
//Generate points 2 steps inside the bounding box.
|
||||||
vtkm::Bounds interiorBounds = bound;
|
vtkm::Bounds interiorBounds = bound;
|
||||||
interiorBounds.X.Min += 2 * stepSize;
|
interiorBounds.X.Min += 2 * stepSize;
|
||||||
@ -467,12 +465,10 @@ void TestEvaluators()
|
|||||||
interiorBounds.X.Max -= 2 * stepSize;
|
interiorBounds.X.Max -= 2 * stepSize;
|
||||||
interiorBounds.Y.Max -= 2 * stepSize;
|
interiorBounds.Y.Max -= 2 * stepSize;
|
||||||
interiorBounds.Z.Max -= 2 * stepSize;
|
interiorBounds.Z.Max -= 2 * stepSize;
|
||||||
for (int k = 0; k < 38; k++)
|
|
||||||
{
|
GenerateRandomParticles(pointIns, 38, interiorBounds);
|
||||||
auto p = RandomPoint(interiorBounds);
|
for (auto& p : pointIns)
|
||||||
pointIns.push_back(vtkm::Particle(p, k));
|
stepResult.push_back(p.Pos + vec * stepSize);
|
||||||
stepResult.push_back(p + vec * stepSize);
|
|
||||||
}
|
|
||||||
|
|
||||||
vtkm::Range xRange, yRange, zRange;
|
vtkm::Range xRange, yRange, zRange;
|
||||||
|
|
||||||
@ -496,8 +492,7 @@ void TestEvaluators()
|
|||||||
// All velocities are in the +ve direction.
|
// All velocities are in the +ve direction.
|
||||||
|
|
||||||
std::vector<vtkm::Particle> boundaryPoints;
|
std::vector<vtkm::Particle> boundaryPoints;
|
||||||
for (int k = 0; k < 10; k++)
|
GenerateRandomParticles(boundaryPoints, 10, forBoundary, 919);
|
||||||
boundaryPoints.push_back(vtkm::Particle(RandomPoint(forBoundary), k));
|
|
||||||
|
|
||||||
for (auto& ds : dataSets)
|
for (auto& ds : dataSets)
|
||||||
{
|
{
|
||||||
@ -506,7 +501,6 @@ void TestEvaluators()
|
|||||||
|
|
||||||
RK4Type rk4(gridEval, stepSize);
|
RK4Type rk4(gridEval, stepSize);
|
||||||
ValidateIntegrator(rk4, pointIns, stepResult, "constant vector RK4");
|
ValidateIntegrator(rk4, pointIns, stepResult, "constant vector RK4");
|
||||||
|
|
||||||
ValidateIntegratorForBoundary(bound, rk4, boundaryPoints, "constant vector RK4");
|
ValidateIntegratorForBoundary(bound, rk4, boundaryPoints, "constant vector RK4");
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
@ -588,9 +582,7 @@ void TestParticleWorkletsWithDataSetTypes()
|
|||||||
|
|
||||||
//Generate three random points.
|
//Generate three random points.
|
||||||
std::vector<vtkm::Particle> pts;
|
std::vector<vtkm::Particle> pts;
|
||||||
pts.push_back(vtkm::Particle(RandomPoint(bound), 0));
|
GenerateRandomParticles(pts, 3, bound, 111);
|
||||||
pts.push_back(vtkm::Particle(RandomPoint(bound), 1));
|
|
||||||
pts.push_back(vtkm::Particle(RandomPoint(bound), 2));
|
|
||||||
std::vector<vtkm::Particle> pts2 = pts;
|
std::vector<vtkm::Particle> pts2 = pts;
|
||||||
|
|
||||||
vtkm::Id nSeeds = static_cast<vtkm::Id>(pts.size());
|
vtkm::Id nSeeds = static_cast<vtkm::Id>(pts.size());
|
||||||
|
Loading…
Reference in New Issue
Block a user