Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-16 17:22:55 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Adding WarpX particle advection streamlines

Browse Source
This commit is contained in:
Abhishek Yenpure 2021-09-01 17:22:51 -07:00
parent e34eafa09b
commit f0d267e813
8 changed files with 83 additions and 79 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

108
vtkm/Particle.h
View File

@ -113,14 +113,6 @@ public:
// troublesome with CUDA __host__ __device__ markup.
}
VTKM_EXEC_CONT
vtkm::Vec3f Next(const vtkm::VecVariable<vtkm::Vec3f, 2>& vectors,
const vtkm::FloatDefault& length)
{
VTKM_ASSERT(vectors.GetNumberOfComponents() > 0);
return this->Pos + length * vectors[0];
}
VTKM_EXEC_CONT
vtkm::Vec3f Velocity(const vtkm::VecVariable<vtkm::Vec3f, 2>& vectors,
const vtkm::FloatDefault& vtkmNotUsed(length))
@ -131,6 +123,13 @@ public:
return vectors[0];
}
VTKM_EXEC_CONT
vtkm::Vec3f GetEvaluationPosition(const vtkm::FloatDefault& deltaT) const
{
(void)deltaT; // unused for a general particle advection case
return this->Pos;
}
inline VTKM_CONT friend std::ostream& operator<<(std::ostream& out, const vtkm::Particle& p)
{
out << "v(" << p.Time << ") = " << p.Pos << ", ID: " << p.ID << ", NumSteps: " << p.NumSteps
@ -145,22 +144,22 @@ public:
vtkm::FloatDefault Time = 0;
};
class Electron
class ChargedParticle
{
public:
VTKM_EXEC_CONT
Electron() {}
ChargedParticle() {}
VTKM_EXEC_CONT
Electron(const vtkm::Vec3f& position,
const vtkm::Id& id,
const vtkm::FloatDefault& mass,
const vtkm::FloatDefault& charge,
const vtkm::FloatDefault& weighting,
const vtkm::Vec3f& momentum,
const vtkm::Id& numSteps = 0,
const vtkm::ParticleStatus& status = vtkm::ParticleStatus(),
const vtkm::FloatDefault& time = 0)
ChargedParticle(const vtkm::Vec3f& position,
const vtkm::Id& id,
const vtkm::FloatDefault& mass,
const vtkm::FloatDefault& charge,
const vtkm::FloatDefault& weighting,
const vtkm::Vec3f& momentum,
const vtkm::Id& numSteps = 0,
const vtkm::ParticleStatus& status = vtkm::ParticleStatus(),
const vtkm::FloatDefault& time = 0)
: Pos(position)
, ID(id)
, NumSteps(numSteps)
@ -173,19 +172,17 @@ public:
{
}
VTKM_EXEC_CONT
vtkm::FloatDefault Beta(vtkm::Vec3f momentum) const
{
return static_cast<vtkm::FloatDefault>(vtkm::Magnitude(momentum / this->Mass) /
vtkm::Pow(SPEED_OF_LIGHT, 2));
}
VTKM_EXEC_CONT
vtkm::Vec3f Next(const vtkm::VecVariable<vtkm::Vec3f, 2>& vectors,
const vtkm::FloatDefault& length)
vtkm::FloatDefault Gamma(vtkm::Vec3f momentum, bool reciprocal = false) const
{
// TODO: implement Lorentz force calculation
return this->Pos + length * this->Velocity(vectors, length);
constexpr vtkm::FloatDefault c2 = SPEED_OF_LIGHT * SPEED_OF_LIGHT;
const auto fMom2 = vtkm::MagnitudeSquared(momentum);
const auto m2 = this->Mass * this->Mass;
const auto m2_c2_reci = 1.0 / (m2 * c2);
if (reciprocal)
return static_cast<vtkm::FloatDefault>(vtkm::RSqrt(1.0 + fMom2 * m2_c2_reci));
else
return static_cast<vtkm::FloatDefault>(vtkm::Sqrt(1.0 + fMom2 * m2_c2_reci));
}
VTKM_EXEC_CONT
@ -197,35 +194,40 @@ public:
// Suppress unused warning
(void)this->Weighting;
vtkm::Vec3f velocity;
// Particle has a charge and a mass
// Velocity updated using Lorentz force
// Return velocity of the particle
vtkm::Vec3f eField = vectors[0];
vtkm::Vec3f bField = vectors[1];
const vtkm::FloatDefault QoM = this->Charge / this->Mass;
const vtkm::FloatDefault half = 0.5;
const vtkm::Vec3f mom_ = this->Momentum + (half * this->Charge * eField * length);
const vtkm::Vec3f mom_minus = this->Momentum + (0.5 * this->Charge * eField * length);
//TODO : Calculate Gamma
vtkm::Vec3f gamma_reci = vtkm::Sqrt(1 - this->Beta(mom_));
// gamma(mom_, mass) -> needs velocity calculation
const vtkm::Vec3f t = half * QoM * bField * gamma_reci * length;
const vtkm::Vec3f s = 2.0f * t * (1 / 1 + vtkm::Magnitude(t));
const vtkm::Vec3f mom_pr = mom_ + vtkm::Cross(mom_, t);
const vtkm::Vec3f mom_pl = mom_ + vtkm::Cross(mom_pr, s);
const vtkm::Vec3f mom_new = mom_pl + 0.5 * this->Charge * eField * length;
//TODO : this is a const method, figure a better way to update momentum
// Get reciprocal of Gamma
vtkm::Vec3f gamma_reci = this->Gamma(mom_minus, true);
const vtkm::Vec3f t = 0.5 * QoM * length * bField * gamma_reci;
const vtkm::Vec3f s = 2.0f * t * (1.0 / (1.0 + vtkm::Magnitude(t)));
const vtkm::Vec3f mom_prime = mom_minus + vtkm::Cross(mom_minus, t);
const vtkm::Vec3f mom_plus = mom_minus + vtkm::Cross(mom_prime, s);
const vtkm::Vec3f mom_new = mom_plus + 0.5 * this->Charge * eField * length;
//TODO : Sould this be a const method?
// If yes, need a better way to update momentum
this->Momentum = mom_new;
velocity = mom_new / this->Mass;
// momentum = velocity * mass * gamma;
// --> velocity = momentum / (mass * gamma)
// --> velocity = ( momentum / mass ) * gamma_reci
vtkm::Vec3f velocity = (mom_new / this->Mass) * this->Gamma(mom_new, true);
return velocity;
}
VTKM_EXEC_CONT
vtkm::Vec3f GetEvaluationPosition(const vtkm::FloatDefault& deltaT) const
{
// Translation is in -ve Z direction,
// this needs to be a parameter.
auto translation = this->NumSteps * deltaT * SPEED_OF_LIGHT * vtkm::Vec3f{ 0., 0., -1.0 };
return this->Pos + translation;
}
vtkm::Vec3f Pos;
vtkm::Id ID = -1;
vtkm::Id NumSteps = 0;
@ -240,7 +242,7 @@ private:
constexpr static vtkm::FloatDefault SPEED_OF_LIGHT =
static_cast<vtkm::FloatDefault>(2.99792458e8);
friend struct mangled_diy_namespace::Serialization<vtkm::Electron>;
friend struct mangled_diy_namespace::Serialization<vtkm::ChargedParticle>;
};
} //namespace vtkm
@ -272,10 +274,10 @@ public:
};
template <>
struct Serialization<vtkm::Electron>
struct Serialization<vtkm::ChargedParticle>
{
public:
static VTKM_CONT void save(BinaryBuffer& bb, const vtkm::Electron& e)
static VTKM_CONT void save(BinaryBuffer& bb, const vtkm::ChargedParticle& e)
{
vtkmdiy::save(bb, e.Pos);
vtkmdiy::save(bb, e.ID);
@ -288,7 +290,7 @@ public:
vtkmdiy::save(bb, e.Momentum);
}
static VTKM_CONT void load(BinaryBuffer& bb, vtkm::Electron& e)
static VTKM_CONT void load(BinaryBuffer& bb, vtkm::ChargedParticle& e)
{
vtkmdiy::load(bb, e.Pos);
vtkmdiy::load(bb, e.ID);

3
vtkm/filter/testing/UnitTestStreamlineFilter.cxx
View File

@ -131,7 +131,7 @@ void TestPathline()
filt.SetNextTime(t1);
filt.SetNextDataSet(ds2);
output = filt.Execute(ds1);
numExpectedPoints = 63;
numExpectedPoints = 33;
}
else
{
@ -149,6 +149,7 @@ void TestPathline()
//Validate the result is correct.
vtkm::cont::CoordinateSystem coords = output.GetCoordinateSystem();
VTKM_TEST_ASSERT(coords.GetNumberOfPoints() == numExpectedPoints,
"Wrong number of coordinates");

2
vtkm/worklet/particleadvection/EulerIntegrator.h
View File

@ -36,7 +36,7 @@ public:
vtkm::Vec3f& velocity) const
{
auto time = particle.Time;
auto inpos = particle.Pos;
auto inpos = particle.GetEvaluationPosition(stepLength);
vtkm::VecVariable<vtkm::Vec3f, 2> vectors;
GridEvaluatorStatus status = this->Evaluator.Evaluate(inpos, time, vectors);
if (status.CheckOk())

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

@ -62,9 +62,9 @@ public:
}
template <typename Point>
VTKM_EXEC bool IsWithinSpatialBoundary(const Point point) const
VTKM_EXEC bool IsWithinSpatialBoundary(const Point& point) const
{
vtkm::Id cellId;
vtkm::Id cellId = -1;
Point parametric;
this->Locator.FindCell(point, cellId, parametric);

9
vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h
View File

@ -53,7 +53,6 @@ public:
const vtkm::Id& maxSteps) const
{
auto particle = integralCurve.GetParticle(idx);
vtkm::FloatDefault time = particle.Time;
bool tookAnySteps = false;
@ -61,16 +60,15 @@ public:
// 1. you could have success AND at temporal boundary.
// 2. could you have success AND at spatial?
// 3. all three?
integralCurve.PreStepUpdate(idx);
do
{
particle = integralCurve.GetParticle(idx);
vtkm::Vec3f outpos;
auto status = integrator.Step(particle, time, outpos);
if (status.CheckOk())
{
integralCurve.StepUpdate(idx, time, outpos);
particle.Pos = outpos;
integralCurve.StepUpdate(idx, particle, time, outpos);
tookAnySteps = true;
}
@ -81,8 +79,7 @@ public:
status = integrator.SmallStep(particle, time, outpos);
if (status.CheckOk())
{
integralCurve.StepUpdate(idx, time, outpos);
particle.Pos = outpos;
integralCurve.StepUpdate(idx, particle, time, outpos);
tookAnySteps = true;
}
}

30
vtkm/worklet/particleadvection/Particles.h
View File

@ -52,13 +52,16 @@ public:
void PreStepUpdate(const vtkm::Id& vtkmNotUsed(idx)) {}
VTKM_EXEC
void StepUpdate(const vtkm::Id& idx, vtkm::FloatDefault time, const vtkm::Vec3f& pt)
void StepUpdate(const vtkm::Id& idx,
const ParticleType& particle,
vtkm::FloatDefault time,
const vtkm::Vec3f& pt)
{
ParticleType p = this->GetParticle(idx);
p.Pos = pt;
p.Time = time;
p.NumSteps++;
this->Particles.Set(idx, p);
ParticleType newParticle(particle);
newParticle.Pos = pt;
newParticle.Time = time;
newParticle.NumSteps++;
this->Particles.Set(idx, newParticle);
}
VTKM_EXEC
@ -66,7 +69,7 @@ public:
const vtkm::worklet::particleadvection::IntegratorStatus& status,
vtkm::Id maxSteps)
{
ParticleType p = this->GetParticle(idx);
ParticleType p(this->GetParticle(idx));
if (p.NumSteps == maxSteps)
p.Status.SetTerminate();
@ -85,7 +88,7 @@ public:
VTKM_EXEC
bool CanContinue(const vtkm::Id& idx)
{
ParticleType p = this->GetParticle(idx);
ParticleType p(this->GetParticle(idx));
return (p.Status.CheckOk() && !p.Status.CheckTerminate() && !p.Status.CheckSpatialBounds() &&
!p.Status.CheckTemporalBounds() && !p.Status.CheckInGhostCell());
@ -94,7 +97,7 @@ public:
VTKM_EXEC
void UpdateTookSteps(const vtkm::Id& idx, bool val)
{
ParticleType p = this->GetParticle(idx);
ParticleType p(this->GetParticle(idx));
if (val)
p.Status.SetTookAnySteps();
else
@ -179,13 +182,14 @@ public:
}
VTKM_EXEC
void StepUpdate(const vtkm::Id& idx, vtkm::FloatDefault time, const vtkm::Vec3f& pt)
void StepUpdate(const vtkm::Id& idx,
const ParticleType& particle,
vtkm::FloatDefault time,
const vtkm::Vec3f& pt)
{
this->ParticleExecutionObject<ParticleType>::StepUpdate(idx, time, pt);
this->ParticleExecutionObject<ParticleType>::StepUpdate(idx, particle, time, pt);
//local step count.
vtkm::Id stepCount = this->StepCount.Get(idx);
vtkm::Id loc = idx * Length + stepCount;
this->History.Set(loc, pt);
this->ValidPoint.Set(loc, 1);

2
vtkm/worklet/particleadvection/RK4Integrator.h
View File

@ -36,7 +36,7 @@ public:
vtkm::Vec3f& velocity) const
{
auto time = particle.Time;
auto inpos = particle.Pos;
auto inpos = particle.GetEvaluationPosition(stepLength);
vtkm::FloatDefault boundary = this->Evaluator.GetTemporalBoundary(static_cast<vtkm::Id>(1));
if ((time + stepLength + vtkm::Epsilon<vtkm::FloatDefault>() - boundary) > 0.0)
stepLength = boundary - time;

4
vtkm/worklet/particleadvection/Stepper.h
View File

@ -74,7 +74,7 @@ public:
}
template <typename Particle>
VTKM_EXEC IntegratorStatus SmallStep(Particle particle,
VTKM_EXEC IntegratorStatus SmallStep(Particle& particle,
vtkm::FloatDefault& time,
vtkm::Vec3f& outpos) const
{
@ -88,7 +88,7 @@ public:
//The binary search will be between {0, this->DeltaT}
vtkm::FloatDefault stepRange[2] = { 0, this->DeltaT };
vtkm::Vec3f currPos(particle.Pos);
vtkm::Vec3f currPos(particle.GetEvaluationPosition(this->DeltaT));
vtkm::Vec3f currVelocity(0, 0, 0);
vtkm::VecVariable<vtkm::Vec3f, 2> currValue, tmp;
auto evalStatus = this->Evaluator.Evaluate(currPos, particle.Time, currValue);