diff --git a/vtkm/exec/CellLocatorRectilinearGrid.h b/vtkm/exec/CellLocatorRectilinearGrid.h index a0eb5803f..14bbb4aa5 100644 --- a/vtkm/exec/CellLocatorRectilinearGrid.h +++ b/vtkm/exec/CellLocatorRectilinearGrid.h @@ -117,6 +117,7 @@ public: if (point[dim] == MaxPoint[dim]) { logicalCell[dim] = this->PointDimensions[dim] - 2; + parametric[dim] = static_cast(1); continue; } diff --git a/vtkm/worklet/particleadvection/GridEvaluators.h b/vtkm/worklet/particleadvection/GridEvaluators.h index af8347b42..d60787c1f 100644 --- a/vtkm/worklet/particleadvection/GridEvaluators.h +++ b/vtkm/worklet/particleadvection/GridEvaluators.h @@ -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(); diff --git a/vtkm/worklet/particleadvection/Integrators.h b/vtkm/worklet/particleadvection/Integrators.h index f7bc01f23..f5bb22336 100644 --- a/vtkm/worklet/particleadvection/Integrators.h +++ b/vtkm/worklet/particleadvection/Integrators.h @@ -13,7 +13,6 @@ #ifndef vtk_m_worklet_particleadvection_Integrators_h #define vtk_m_worklet_particleadvection_Integrators_h -#include #include #include @@ -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(); + const vtkm::FloatDefault eps = vtkm::Epsilon() * 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() - 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(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(6); return IntegratorStatus(true, false, evalStatus.CheckTemporalBounds()); } }; diff --git a/vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h b/vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h index c6221f3ea..43db80967 100644 --- a/vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h +++ b/vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h @@ -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<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< #include -#include - -//timers -#include -#include - - +#include namespace { @@ -98,20 +92,29 @@ vtkm::FloatDefault vecData[125 * 3] = { }; } -vtkm::Vec3f RandomPoint(const vtkm::Bounds& bounds) +void GenerateRandomParticles(std::vector& points, + std::size_t N, + vtkm::Bounds& bounds, + std::size_t seed = 314) { - vtkm::FloatDefault rx = - static_cast(rand()) / static_cast(RAND_MAX); - vtkm::FloatDefault ry = - static_cast(rand()) / static_cast(RAND_MAX); - vtkm::FloatDefault rz = - static_cast(rand()) / static_cast(RAND_MAX); + std::random_device device; + std::default_random_engine generator(static_cast(seed)); + vtkm::FloatDefault zero(0), one(1); + std::uniform_real_distribution distribution(zero, one); - vtkm::Vec3f p; - p[0] = static_cast(bounds.X.Min + rx * bounds.X.Length()); - p[1] = static_cast(bounds.Y.Min + ry * bounds.Y.Length()); - p[2] = static_cast(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(bounds.X.Min + rx * bounds.X.Length()); + p[1] = static_cast(bounds.Y.Min + ry * bounds.Y.Length()); + p[2] = static_cast(bounds.Z.Min + rz * bounds.Z.Length()); + points.push_back(vtkm::Particle(p, static_cast(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(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 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 pointIns; std::vector 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 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 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 pts2 = pts; vtkm::Id nSeeds = static_cast(pts.size());