mirror of
https://gitlab.kitware.com/vtk/vtk-m
synced 2024-09-20 02:55:47 +00:00
Particle Advection Enhancements
-Added support for short steps with a generic way to support for al types of integrators -Added a base class for all integrators to leverage the generic way to handle taking steps and push particle out of spatial domain
This commit is contained in:
Yenpure
parent
12d8705796
commit
a3938245d2
@ -78,27 +78,47 @@ public:
|
||||
VTKM_EXEC
|
||||
virtual FieldType GetEscapeStepLength(const vtkm::Vec<FieldType, 3>& inpos,
|
||||
const PortalType& field,
|
||||
FieldType stepLength) const = 0;
|
||||
vtkm::Vec<FieldType, 3>& velocity,
|
||||
FieldType stepLength) const
|
||||
{
|
||||
this->CheckStep(inpos, field, stepLength, velocity);
|
||||
FieldType magnitude = vtkm::Magnitude(velocity);
|
||||
vtkm::Vec<FieldType, 3> dir = velocity / magnitude;
|
||||
FieldType xbound = static_cast<FieldType>(dir[0] > 0 ? Bounds.X.Max : Bounds.X.Min);
|
||||
FieldType ybound = static_cast<FieldType>(dir[1] > 0 ? Bounds.Y.Max : Bounds.Y.Min);
|
||||
FieldType zbound = static_cast<FieldType>(dir[2] > 0 ? Bounds.Z.Max : Bounds.Z.Min);
|
||||
/*Add a fraction just push the particle beyond the bounds*/
|
||||
FieldType hx = (std::abs(xbound - inpos[0]) + this->Tolerance) / std::abs(velocity[0]);
|
||||
FieldType hy = (std::abs(ybound - inpos[1]) + this->Tolerance) / std::abs(velocity[1]);
|
||||
FieldType hz = (std::abs(zbound - inpos[2]) + this->Tolerance) / std::abs(velocity[2]);
|
||||
stepLength = std::min(hx, std::min(hy, hz));
|
||||
return stepLength;
|
||||
}
|
||||
|
||||
VTKM_EXEC
|
||||
ParticleStatus PushOutOfDomain(vtkm::Vec<FieldType, 3> inpos,
|
||||
const PortalType& field,
|
||||
vtkm::Id numSteps,
|
||||
vtkm::Vec<FieldType, 3>& outpos) const
|
||||
{
|
||||
numSteps = (numSteps == 0) ? 1 : numSteps;
|
||||
FieldType totalTime = numSteps * this->StepLength;
|
||||
FieldType timeFraction = totalTime * this->Tolerance;
|
||||
FieldType stepLength = this->StepLength / 2;
|
||||
vtkm::Vec<FieldType, 3> velocity;
|
||||
vtkm::Id numSteps = 0;
|
||||
do
|
||||
if (this->ShortStepsSupported)
|
||||
{
|
||||
ParticleStatus status = this->CheckStep(inpos, field, stepLength, velocity);
|
||||
if (status == ParticleStatus::STATUS_OK)
|
||||
do
|
||||
{
|
||||
inpos = inpos + stepLength * velocity;
|
||||
++numSteps;
|
||||
}
|
||||
stepLength = stepLength / 2;
|
||||
} while (numSteps < 10);
|
||||
stepLength = GetEscapeStepLength(inpos, field, stepLength);
|
||||
ParticleStatus status = this->CheckStep(inpos, field, stepLength, velocity);
|
||||
if (status == ParticleStatus::STATUS_OK)
|
||||
{
|
||||
inpos = inpos + stepLength * velocity;
|
||||
}
|
||||
stepLength = stepLength / 2;
|
||||
} while (stepLength > timeFraction);
|
||||
}
|
||||
stepLength = GetEscapeStepLength(inpos, field, velocity, stepLength);
|
||||
outpos = inpos + stepLength * velocity;
|
||||
return ParticleStatus::EXITED_SPATIAL_BOUNDARY;
|
||||
}
|
||||
@ -106,6 +126,9 @@ public:
|
||||
protected:
|
||||
FieldEvaluateType Evaluator;
|
||||
FieldType StepLength;
|
||||
FieldType Tolerance = static_cast<FieldType>(1e-6);
|
||||
bool ShortStepsSupported = false;
|
||||
vtkm::Bounds Bounds;
|
||||
};
|
||||
|
||||
template <typename FieldEvaluateType, typename FieldType, typename PortalType>
|
||||
@ -116,13 +139,15 @@ public:
|
||||
RK4Integrator()
|
||||
: Integrator<FieldEvaluateType, FieldType, PortalType>()
|
||||
{
|
||||
this->ShortStepsSupported = true;
|
||||
}
|
||||
|
||||
VTKM_EXEC_CONT
|
||||
RK4Integrator(const FieldEvaluateType& field, FieldType stepLength, vtkm::cont::DataSet& dataset)
|
||||
: Integrator<FieldEvaluateType, FieldType, PortalType>(field, stepLength)
|
||||
{
|
||||
bounds = dataset.GetCoordinateSystem(0).GetBounds();
|
||||
this->ShortStepsSupported = true;
|
||||
this->Bounds = dataset.GetCoordinateSystem(0).GetBounds();
|
||||
}
|
||||
|
||||
VTKM_EXEC
|
||||
@ -131,7 +156,7 @@ public:
|
||||
FieldType stepLength,
|
||||
vtkm::Vec<FieldType, 3>& velocity) const VTKM_OVERRIDE
|
||||
{
|
||||
if (!bounds.Contains(inpos))
|
||||
if (!this->Bounds.Contains(inpos))
|
||||
{
|
||||
return ParticleStatus::EXITED_SPATIAL_BOUNDARY;
|
||||
}
|
||||
@ -150,163 +175,33 @@ public:
|
||||
return ParticleStatus::AT_SPATIAL_BOUNDARY;
|
||||
}
|
||||
}
|
||||
|
||||
VTKM_EXEC
|
||||
virtual FieldType GetEscapeStepLength(const vtkm::Vec<FieldType, 3>& inpos,
|
||||
const PortalType& field,
|
||||
FieldType stepLength) const VTKM_OVERRIDE
|
||||
{
|
||||
vtkm::Vec<FieldType, 3> velocity;
|
||||
this->CheckStep(inpos, field, stepLength, velocity);
|
||||
FieldType magnitude = vtkm::Magnitude(velocity);
|
||||
vtkm::Vec<FieldType, 3> dir = velocity / magnitude;
|
||||
FieldType xbound = static_cast<FieldType>(dir[0] > 0 ? bounds.X.Max : bounds.X.Min);
|
||||
FieldType ybound = static_cast<FieldType>(dir[1] > 0 ? bounds.Y.Max : bounds.Y.Min);
|
||||
FieldType zbound = static_cast<FieldType>(dir[2] > 0 ? bounds.Z.Max : bounds.Z.Min);
|
||||
FieldType hx = std::abs(xbound - inpos[0]) / std::abs(velocity[0]);
|
||||
FieldType hy = std::abs(ybound - inpos[1]) / std::abs(velocity[1]);
|
||||
FieldType hz = std::abs(zbound - inpos[2]) / std::abs(velocity[2]);
|
||||
stepLength = std::min(hx, std::min(hy, hz));
|
||||
stepLength += stepLength / 100.0f;
|
||||
return stepLength;
|
||||
}
|
||||
|
||||
|
||||
private:
|
||||
vtkm::Bounds bounds;
|
||||
};
|
||||
|
||||
/*template<typename FieldEvaluateType, typename FieldType>
|
||||
class RK4Integrator
|
||||
template <typename FieldEvaluateType, typename FieldType, typename PortalType>
|
||||
class EulerIntegrator : public Integrator<FieldEvaluateType, FieldType, PortalType>
|
||||
{
|
||||
public:
|
||||
VTKM_EXEC_CONT
|
||||
RK4Integrator()
|
||||
: h(0)
|
||||
EulerIntegrator(const FieldEvaluateType& evaluator, FieldType field, vtkm::cont::DataSet& dataset)
|
||||
: Integrator<FieldEvaluateType, FieldType, PortalType>(evaluator, field)
|
||||
{
|
||||
this->ShortStepsSupported = false;
|
||||
this->Bounds = dataset.GetCoordinateSystem(0).GetBounds();
|
||||
}
|
||||
|
||||
VTKM_EXEC_CONT
|
||||
RK4Integrator(const FieldEvaluateType& field, FieldType _h, vtkm::cont::DataSet& dataset)
|
||||
: f(field)
|
||||
, h(_h)
|
||||
VTKM_EXEC
|
||||
virtual ParticleStatus CheckStep(const vtkm::Vec<FieldType, 3>& inpos,
|
||||
const PortalType& field,
|
||||
vtkm::Vec<FieldType, 3>& velocity) const
|
||||
{
|
||||
bounds = dataset.GetCoordinateSystem(0).GetBounds();
|
||||
}
|
||||
|
||||
template<typename PortalType>
|
||||
VTKM_EXEC ParticleStatus Step(const vtkm::Vec<FieldType, 3>& pos,
|
||||
const PortalType& field,
|
||||
vtkm::Vec<FieldType, 3>& out) const
|
||||
{
|
||||
if(!bounds.Contains(pos))
|
||||
{
|
||||
out[0] = pos[0]; out[1] = pos[1]; out[2] = pos[2];
|
||||
return ParticleStatus::EXITED_SPATIAL_BOUNDARY;
|
||||
}
|
||||
vtkm::Vec<FieldType, 3> k1, k2, k3, k4;
|
||||
vtkm::Id numShortSteps = 0;
|
||||
bool shortSteps = false;
|
||||
FieldType step_h = h;
|
||||
FieldType step_h_2 = step_h / 2;
|
||||
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);
|
||||
bool validPos = this->Evaluator.Evaluate(inpos, field, velocity);
|
||||
if (validPos)
|
||||
return ParticleStatus::STATUS_OK;
|
||||
}
|
||||
else
|
||||
{
|
||||
shortSteps = true;
|
||||
}
|
||||
*Take short steps to minimize advection error*
|
||||
while (shortSteps)
|
||||
{
|
||||
*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) &&
|
||||
f.Evaluate(pos + step_h_2 * k2, field, k3) && f.Evaluate(pos + step_h * k3, field, k4))
|
||||
{
|
||||
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*
|
||||
if (numShortSteps == 2)
|
||||
{
|
||||
step_h /= 2;
|
||||
step_h_2 = step_h / 2;
|
||||
*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<FieldType, 3> vel = (k1 + 2 * k2 + 2 * k3 + k4) / 6.0f;
|
||||
*Get the direction of the particle*
|
||||
FieldType magnitude = vtkm::Magnitude(vel);
|
||||
vtkm::Vec<FieldType, 3> dir = vel / magnitude;
|
||||
*Determine the bounds for the particle*
|
||||
FieldType xbound = static_cast<FieldType>(dir[0] > 0 ? bounds.X.Max : bounds.X.Min);
|
||||
FieldType ybound = static_cast<FieldType>(dir[1] > 0 ? bounds.Y.Max : bounds.Y.Min);
|
||||
FieldType zbound = static_cast<FieldType>(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]);
|
||||
FieldType hesc = std::min(hx, std::min(hy, hz));
|
||||
hesc += hesc / 100.0f;
|
||||
out += hesc * vel;
|
||||
shortSteps = false;
|
||||
break;
|
||||
}
|
||||
}
|
||||
return ParticleStatus::EXITED_SPATIAL_BOUNDARY;
|
||||
return ParticleStatus::AT_SPATIAL_BOUNDARY;
|
||||
}
|
||||
|
||||
FieldEvaluateType f;
|
||||
FieldType h;
|
||||
vtkm::Bounds bounds;
|
||||
};*/ //RK4Integrator
|
||||
|
||||
/*template<typename FieldEvaluateType, typename FieldType>
|
||||
class EulerIntegrator
|
||||
{
|
||||
public:
|
||||
VTKM_EXEC_CONT
|
||||
EulerIntegrator(const FieldEvaluateType& field, FieldType _h)
|
||||
: f(field)
|
||||
, h(_h)
|
||||
{
|
||||
}
|
||||
|
||||
template<typename PortalType>
|
||||
VTKM_EXEC ParticleStatus Step(const vtkm::Vec<FieldType, 3>& pos,
|
||||
const PortalType& field,
|
||||
vtkm::Vec<FieldType, 3>& out) const
|
||||
{
|
||||
vtkm::Vec<FieldType, 3> vCur;
|
||||
if (f.Evaluate(pos, field, vCur))
|
||||
{
|
||||
out = pos + h * vCur;
|
||||
return ParticleStatus::STATUS::OK;
|
||||
}
|
||||
return ParticleStatus::STATUS::EXITED_SPATIAL_BOUNDARY;
|
||||
}
|
||||
|
||||
FieldEvaluateType f;
|
||||
FieldType h;
|
||||
};*/ //EulerIntegrator
|
||||
}; //EulerIntegrator
|
||||
|
||||
} //namespace particleadvection
|
||||
} //namespace worklet
|
||||
|
||||
@ -66,7 +66,8 @@ public:
|
||||
}
|
||||
if (status == ParticleStatus::AT_SPATIAL_BOUNDARY)
|
||||
{
|
||||
status = integrator.PushOutOfDomain(inpos, field, outpos);
|
||||
vtkm::Id numSteps = ic.GetStep(idx);
|
||||
status = integrator.PushOutOfDomain(inpos, field, numSteps, outpos);
|
||||
}
|
||||
if (status == ParticleStatus::EXITED_SPATIAL_BOUNDARY)
|
||||
{
|
||||
|
||||
Loading…
Reference in New Issue
Block a user