From 61c61ee93bce85c63dbb21ffaa70688bdf7d2f51 Mon Sep 17 00:00:00 2001 From: Yenpure Date: Mon, 31 Jul 2017 15:46:13 -0400 Subject: [PATCH] Adding a base class for Integrators --- .../particle_advection/ParticleAdvection.cxx | 5 +- vtkm/worklet/particleadvection/Integrators.h | 134 ++++++++++++++---- .../ParticleAdvectionWorklets.h | 19 ++- vtkm/worklet/particleadvection/Particles.h | 8 +- .../testing/UnitTestParticleAdvection.cxx | 5 +- 5 files changed, 128 insertions(+), 43 deletions(-) diff --git a/examples/particle_advection/ParticleAdvection.cxx b/examples/particle_advection/ParticleAdvection.cxx index 78659aecb..6f4c0cfe0 100644 --- a/examples/particle_advection/ParticleAdvection.cxx +++ b/examples/particle_advection/ParticleAdvection.cxx @@ -76,7 +76,10 @@ void RunTest(const std::string& fname, typedef vtkm::worklet::particleadvection::UniformGridEvaluate RGEvalType; - typedef vtkm::worklet::particleadvection::RK4Integrator RK4RGType; + typedef vtkm::worklet::particleadvection::RK4Integrator + RK4RGType; RGEvalType eval(ds); RK4RGType rk4(eval, stepSize, ds); diff --git a/vtkm/worklet/particleadvection/Integrators.h b/vtkm/worklet/particleadvection/Integrators.h index c331b8381..1fbce7434 100644 --- a/vtkm/worklet/particleadvection/Integrators.h +++ b/vtkm/worklet/particleadvection/Integrators.h @@ -40,7 +40,6 @@ protected: VTKM_EXEC_CONT Integrator() : StepLength(0) - , HighAccuracy(false) { } @@ -48,33 +47,108 @@ protected: Integrator(const FieldEvaluateType& evaluator, FieldType stepLength) : Evaluator(evaluator) , StepLength(stepLength) - , HighAccuracy(false) { } +public: VTKM_EXEC - virtual ParticleStatus Step(const vtkm::Vec& inpos, - const PortalType& vectorField, - vtkm::Vec& outpos) const; + ParticleStatus Step(const vtkm::Vec& inpos, + const PortalType& field, + vtkm::Vec& outpos) const + { + vtkm::Vec velocity; + ParticleStatus status = this->CheckStep(inpos, field, this->StepLength, velocity); + if (status == ParticleStatus::STATUS_OK) + { + outpos = inpos + this->StepLength * velocity; + } + else + { + outpos = inpos; + } + return status; + } VTKM_EXEC - virtual bool IsStepValid(const vtkm::Vec& inpos, - const PortalType& vectorField, - vtkm::Vec& velocity) const; + virtual ParticleStatus CheckStep(const vtkm::Vec& inpos, + const PortalType& field, + FieldType stepLength, + vtkm::Vec& velocity) const; VTKM_EXEC - virtual ParticleStatus PushOutOfDomain(const vtkm::Vec& inpos, - const PortalType& vectorField, - vtkm::Vec& outpos) const; + ParticleStatus PushOutOfDomain(const vtkm::Vec& inpos, + const PortalType& field, + vtkm::Vec& outpos) const + { + FieldType stepLength = this->StepLength; + vtkm::Vec velocity; + vtkm::Id numSteps = 0; + do + { + stepLength = stepLength / 2; + ParticleStatus status = this->CheckStep(inpos, field, stepLength, velocity); + if (status == ParticleStatus::STATUS_OK) + { + outpos = inpos + stepLength * velocity; + ++numSteps; + } + } while (numSteps < 10); + return ParticleStatus::EXITED_SPATIAL_BOUNDARY; + } protected: FieldEvaluateType Evaluator; FieldType StepLength; - bool HighAccuracy; }; +template +class RK4Integrator : public Integrator +{ +public: + VTKM_EXEC_CONT + RK4Integrator() + : Integrator() + { + } -template + VTKM_EXEC_CONT + RK4Integrator(const FieldEvaluateType& field, FieldType stepLength, vtkm::cont::DataSet& dataset) + : Integrator(field, stepLength) + { + bounds = dataset.GetCoordinateSystem(0).GetBounds(); + } + + VTKM_EXEC + ParticleStatus CheckStep(const vtkm::Vec& inpos, + const PortalType& field, + FieldType stepLength, + vtkm::Vec& velocity) const + { + if (!bounds.Contains(inpos)) + { + return ParticleStatus::EXITED_SPATIAL_BOUNDARY; + } + vtkm::Vec k1, k2, k3, k4; + bool firstOrderValid = this->Evaluator.Evaluate(inpos, field, k1); + bool secondOrderValid = this->Evaluator.Evaluate(inpos + (stepLength / 2) * k1, field, k2); + bool thirdOrderValid = this->Evaluator.Evaluate(inpos + (stepLength / 2) * k2, field, k3); + bool fourthOrderValid = this->Evaluator.Evaluate(inpos + stepLength * k3, field, k4); + velocity = stepLength / 6.0f * (k1 + 2 * k2 + 2 * k3 + k4); + if (firstOrderValid && secondOrderValid && thirdOrderValid && fourthOrderValid) + { + return ParticleStatus::STATUS_OK; + } + else + { + return ParticleStatus::AT_SPATIAL_BOUNDARY; + } + } + +private: + vtkm::Bounds bounds; +}; + +/*template class RK4Integrator { public: @@ -92,16 +166,14 @@ public: bounds = dataset.GetCoordinateSystem(0).GetBounds(); } - template + template VTKM_EXEC ParticleStatus Step(const vtkm::Vec& pos, const PortalType& field, vtkm::Vec& out) const { - if (!bounds.Contains(pos)) + if(!bounds.Contains(pos)) { - out[0] = pos[0]; - out[1] = pos[1]; - out[2] = pos[2]; + out[0] = pos[0]; out[1] = pos[1]; out[2] = pos[2]; return ParticleStatus::EXITED_SPATIAL_BOUNDARY; } vtkm::Vec k1, k2, k3, k4; @@ -112,10 +184,10 @@ public: if (f.Evaluate(pos, field, k1) && f.Evaluate(pos + step_h_2 * k1, field, k2) && f.Evaluate(pos + step_h_2 * k2, field, k3) && f.Evaluate(pos + step_h * k3, field, k4)) { - /* + * * If the particle is inside bounds after taking the steps * return that the step was successful. - */ + * out = pos + step_h / 6.0f * (k1 + 2 * k2 + 2 * k3 + k4); return ParticleStatus::STATUS_OK; } @@ -123,10 +195,10 @@ public: { shortSteps = true; } - /*Take short steps to minimize advection error*/ + *Take short steps to minimize advection error* while (shortSteps) { - /*reduce the step size to half*/ + *reduce the step size to half* step_h /= 2; step_h_2 = step_h / 2; if (f.Evaluate(pos, field, k1) && f.Evaluate(pos + step_h_2 * k1, field, k2) && @@ -135,32 +207,32 @@ public: out = pos + step_h / 6.0f * (k1 + 2 * k2 + 2 * k3 + k4); numShortSteps++; } - /* + * * At this time the particle is still inside the bounds * To Push it at/out of the boundary take an Euler step - */ - /*Check for the function like VisIt*/ + * + *Check for the function like VisIt* if (numShortSteps == 2) { step_h /= 2; step_h_2 = step_h / 2; - /*Calculate the velocity of the particle at current position*/ + *Calculate the velocity of the particle at current position* f.Evaluate(out, field, k1); f.Evaluate(out + step_h_2 * k1, field, k2); f.Evaluate(out + step_h_2 * k2, field, k3); f.Evaluate(out + step_h * k3, field, k4); vtkm::Vec vel = (k1 + 2 * k2 + 2 * k3 + k4) / 6.0f; - /*Get the direction of the particle*/ + *Get the direction of the particle* FieldType magnitude = vtkm::Magnitude(vel); vtkm::Vec dir = vel / magnitude; - /*Determine the bounds for the particle*/ + *Determine the bounds for the particle* FieldType xbound = static_cast(dir[0] > 0 ? bounds.X.Max : bounds.X.Min); FieldType ybound = static_cast(dir[1] > 0 ? bounds.Y.Max : bounds.Y.Min); FieldType zbound = static_cast(dir[2] > 0 ? bounds.Z.Max : bounds.Z.Min); - /* + * * Determine the minimum travel time for the * particle to reach the boundary - */ + * FieldType hx = std::abs(xbound - out[0]) / std::abs(vel[0]); FieldType hy = std::abs(ybound - out[1]) / std::abs(vel[1]); FieldType hz = std::abs(zbound - out[2]) / std::abs(vel[2]); @@ -177,7 +249,7 @@ public: FieldEvaluateType f; FieldType h; vtkm::Bounds bounds; -}; //RK4Integrator +};*/ //RK4Integrator /*template class EulerIntegrator diff --git a/vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h b/vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h index 4619e8bac..49296708a 100644 --- a/vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h +++ b/vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h @@ -53,20 +53,25 @@ public: template VTKM_EXEC void operator()(const vtkm::Id& idx, IntegralCurveType& ic) const { - vtkm::Vec p = ic.GetPos(idx); - vtkm::Vec p2; + vtkm::Vec inpos = ic.GetPos(idx); + vtkm::Vec outpos; while (!ic.Done(idx)) { - ParticleStatus status = integrator.Step(p, field, p2); + ParticleStatus status = integrator.Step(inpos, field, outpos); if (status == ParticleStatus::STATUS_OK) { - ic.TakeStep(idx, p2, status); - p = p2; + ic.TakeStep(idx, outpos, status); + inpos = outpos; } - else if (status == ParticleStatus::EXITED_SPATIAL_BOUNDARY) + if (status == ParticleStatus::AT_SPATIAL_BOUNDARY) { - ic.TakeStep(idx, p2, status); + status = integrator.PushOutOfDomain(inpos, field, outpos); + } + if (status == ParticleStatus::EXITED_SPATIAL_BOUNDARY) + { + ic.TakeStep(idx, outpos, status); + ic.SetExitedSpatialBoundary(idx); } } } diff --git a/vtkm/worklet/particleadvection/Particles.h b/vtkm/worklet/particleadvection/Particles.h index 52050b5cc..3b52889a9 100644 --- a/vtkm/worklet/particleadvection/Particles.h +++ b/vtkm/worklet/particleadvection/Particles.h @@ -36,9 +36,11 @@ enum ParticleStatus { STATUS_OK = 1, TERMINATED = 1 << 1, - EXITED_SPATIAL_BOUNDARY = 1 << 2, - EXITED_TEMPORAL_BOUNDARY = 1 << 3, - STATUS_ERROR = 1 << 4 + AT_SPATIAL_BOUNDARY = 1 << 2, + AT_TEMPORAL_BOUNDARY = 1 << 3, + EXITED_SPATIAL_BOUNDARY = 1 << 4, + EXITED_TEMPORAL_BOUNDARY = 1 << 5, + STATUS_ERROR = 1 << 6 }; template diff --git a/vtkm/worklet/testing/UnitTestParticleAdvection.cxx b/vtkm/worklet/testing/UnitTestParticleAdvection.cxx index c01249d06..4c2be17c3 100644 --- a/vtkm/worklet/testing/UnitTestParticleAdvection.cxx +++ b/vtkm/worklet/testing/UnitTestParticleAdvection.cxx @@ -128,7 +128,10 @@ void TestParticleAdvection() typedef vtkm::worklet::particleadvection::UniformGridEvaluate RGEvalType; - typedef vtkm::worklet::particleadvection::RK4Integrator RK4RGType; + typedef vtkm::worklet::particleadvection::RK4Integrator + RK4RGType; RGEvalType eval(ds); RK4RGType rk4(eval, stepSize, ds);