Shuenhoy/vtk-m2
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Fixing problems with Streamline unit test

Browse Source 
- Fixing issues with temporal bounds to take a step large enough to end
  advection beyond which it is no longer possible
This commit is contained in:
ayenpure 2019-08-01 11:45:28 -06:00
parent d636dea4e3
commit d55b9e3eb5
3 changed files with 17 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -69,21 +69,14 @@ public:
bool IsWithinTemporalBoundary(const vtkm::FloatDefault vtkmNotUsed(time)) const { return true; }
VTKM_EXEC
void GetSpatialBoundary(vtkm::Bounds& bounds) const
{
// Based on the direction of the velocity we need to be able to tell where
// the particle will exit the domain from to actually push it out of domain.
/*boundary[0] = static_cast<ScalarType>(dir[0] > 0 ? this->Bounds.X.Max : this->Bounds.X.Min);
boundary[1] = static_cast<ScalarType>(dir[1] > 0 ? this->Bounds.Y.Max : this->Bounds.Y.Min);
boundary[2] = static_cast<ScalarType>(dir[2] > 0 ? this->Bounds.Z.Max : this->Bounds.Z.Min);*/
bounds = this->Bounds;
}
vtkm::Bounds GetSpatialBoundary() const { return this->Bounds; }
VTKM_EXEC_CONT
void GetTemporalBoundary(vtkm::FloatDefault& boundary) const
vtkm::FloatDefault GetTemporalBoundary(vtkm::Id direction) const
{
// Return the time of the newest time slice
boundary = 0;
return direction > 0 ? vtkm::Infinity<vtkm::FloatDefault>()
: vtkm::NegativeInfinity<vtkm::FloatDefault>();
}
template <typename Point>

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

@ -109,8 +109,8 @@ protected:
ScalarType& time,
vtkm::Vec<ScalarType, 3>& outpos) const override
{
// If without taking the step the particle is out of either spatial
// or temporal boundary, then return the corresponding status.
// If particle is out of either spatial or temporal boundary to begin with,
// then return the corresponding status.
if (!this->Evaluator.IsWithinSpatialBoundary(inpos))
return ParticleStatus::AT_SPATIAL_BOUNDARY;
if (!this->Evaluator.IsWithinTemporalBoundary(time))
@ -139,7 +139,6 @@ protected:
return ParticleStatus::AT_SPATIAL_BOUNDARY;
if (!this->Evaluator.IsWithinTemporalBoundary(time))
return ParticleStatus::AT_TEMPORAL_BOUNDARY;
const vtkm::Float64 eps = 1e-6;
ScalarType optimalLength = static_cast<ScalarType>(0);
vtkm::Id iteration = static_cast<vtkm::Id>(1);
vtkm::Id maxIterations = static_cast<vtkm::Id>(1 << 20);
@ -168,9 +167,10 @@ protected:
// using the higher order evaluator, take a step using that evaluator.
// Take one final step, which should be an Euler step just to push the
// particle out of the domain boundary
vtkm::Bounds bounds;
this->Evaluator.GetSpatialBoundary(bounds);
vtkm::Bounds bounds = this->Evaluator.GetSpatialBoundary();
vtkm::Vec<ScalarType, 3> direction = velocity / vtkm::Magnitude(velocity);
const vtkm::Float64 eps = 1e-6;
ScalarType xStepLength = vtkm::Abs(direction[0] * eps * bounds.X.Length());
ScalarType yStepLength = vtkm::Abs(direction[1] * eps * bounds.Y.Length());
ScalarType zStepLength = vtkm::Abs(direction[2] * eps * bounds.Z.Length());
@ -260,6 +260,10 @@ public:
ScalarType time,
vtkm::Vec<ScalarType, 3>& velocity) const
{
ScalarType boundary = this->Evaluator.GetTemporalBoundary(static_cast<vtkm::Id>(1));
if ((time + stepLength + vtkm::Epsilon<ScalarType>() - boundary) > 0.0)
stepLength = boundary - time;
ScalarType var1 = (stepLength / static_cast<ScalarType>(2));
ScalarType var2 = time + var1;
ScalarType var3 = time + stepLength;
@ -345,11 +349,7 @@ public:
ScalarType time,
vtkm::Vec<ScalarType, 3>& velocity) const
{
bool result = this->Evaluator.Evaluate(inpos, time, velocity);
if (result)
return ParticleStatus::STATUS_OK;
else
return ParticleStatus::EXITED_SPATIAL_BOUNDARY;
return this->Evaluator.Evaluate(inpos, time, velocity);
}
};

12
vtkm/worklet/particleadvection/TemporalGridEvaluators.h
View File

@ -59,18 +59,12 @@ public:
}
VTKM_EXEC
void GetSpatialBoundary(vtkm::Bounds& bounds) const
{
// Based on the direction of the velocity we need to be able to tell where
// the particle will exit the domain from to actually push it out of domain.
return this->EvaluatorTwo.GetSpatialBoundary(bounds);
}
vtkm::Bounds GetSpatialBoundary() const { return this->EvaluatorTwo.GetSpatialBoundary(); }
VTKM_EXEC_CONT
void GetTemporalBoundary(vtkm::FloatDefault& boundary) const
vtkm::FloatDefault GetTemporalBoundary(vtkm::Id direction) const
{
// Return the time of the newest time slice
boundary = TimeTwo;
return direction > 0 ? this->TimeTwo : this->TimeOne;
}
template <typename Point>