Merge topic 'refactor_integrators'

8124fe81f Add missing #include.
6ea07cc1e Missing #include.
024ab1cff Put each ODE integrator in it's own file.

Acked-by: Kitware Robot <kwrobot@kitware.com>
Acked-by: Dave Pugmire <dpugmire@gmail.com>
Acked-by: Kenneth Moreland <kmorel@sandia.gov>
Merge-request: !2213
This commit is contained in:
Nick Thompson 2020-08-18 19:36:49 +00:00 committed by Kitware Robot
commit 88a0cff7cb
20 changed files with 267 additions and 202 deletions

@ -25,8 +25,9 @@
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/particleadvection/Field.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
#include <vtkm/worklet/particleadvection/Particles.h>
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
#include <cstring>
#include <sstream>

@ -17,7 +17,7 @@
#include <vtkm/worklet/ParticleAdvection.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
namespace vtkm
{

@ -16,8 +16,9 @@
#include <vtkm/cont/Invoker.h>
#include <vtkm/worklet/particleadvection/Field.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
#include <vtkm/worklet/particleadvection/Particles.h>
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
#include <vtkm/worklet/LagrangianStructures.h>

@ -14,7 +14,7 @@
#include <vtkm/filter/FilterDataSetWithField.h>
#include <vtkm/worklet/ParticleAdvection.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
namespace vtkm
{

@ -18,8 +18,8 @@
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/cont/ParticleArrayCopy.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/Particles.h>
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
namespace vtkm
{

@ -14,7 +14,7 @@
#include <vtkm/filter/FilterDataSetWithField.h>
#include <vtkm/worklet/ParticleAdvection.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
namespace vtkm
{

@ -16,8 +16,9 @@
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/worklet/particleadvection/Field.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
#include <vtkm/worklet/particleadvection/Particles.h>
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
#include <vtkm/worklet/particleadvection/TemporalGridEvaluators.h>
namespace vtkm

@ -15,7 +15,7 @@
#include <vtkm/worklet/ParticleAdvection.h>
#include <vtkm/worklet/StreamSurface.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
namespace vtkm
{

@ -19,8 +19,8 @@
#include <vtkm/worklet/ParticleAdvection.h>
#include <vtkm/worklet/particleadvection/Field.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/Particles.h>
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
namespace vtkm
{

@ -14,7 +14,8 @@
#include <vtkm/filter/FilterDataSetWithField.h>
#include <vtkm/worklet/ParticleAdvection.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
namespace vtkm
{

@ -17,7 +17,7 @@
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/worklet/particleadvection/Field.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
#include <vtkm/worklet/particleadvection/Particles.h>
namespace vtkm

@ -10,13 +10,15 @@
set(headers
CellInterpolationHelper.h
EulerIntegrator.h
Field.h
GridEvaluators.h
GridEvaluatorStatus.h
Integrators.h
IntegratorBase.h
IntegratorStatus.h
Particles.h
ParticleAdvectionWorklets.h
RK4Integrator.h
TemporalGridEvaluators.h
)

@ -0,0 +1,97 @@
//=============================================================================
//
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
//=============================================================================
#ifndef vtk_m_worklet_particleadvection_EulerIntegrator_h
#define vtk_m_worklet_particleadvection_EulerIntegrator_h
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
namespace vtkm
{
namespace worklet
{
namespace particleadvection
{
template <typename FieldEvaluateType>
class EulerIntegrator : public IntegratorBase
{
public:
EulerIntegrator() = default;
VTKM_CONT
EulerIntegrator(const FieldEvaluateType& evaluator, const vtkm::FloatDefault stepLength)
: IntegratorBase(stepLength)
, Evaluator(evaluator)
{
}
template <typename Device>
class ExecObject : public IntegratorBase::ExecObjectBaseImpl<
vtkm::cont::internal::ExecutionObjectType<FieldEvaluateType, Device>,
typename EulerIntegrator::template ExecObject<Device>>
{
VTKM_IS_DEVICE_ADAPTER_TAG(Device);
using FieldEvaluateExecType =
vtkm::cont::internal::ExecutionObjectType<FieldEvaluateType, Device>;
using Superclass =
IntegratorBase::ExecObjectBaseImpl<FieldEvaluateExecType,
typename EulerIntegrator::template ExecObject<Device>>;
public:
VTKM_EXEC_CONT
ExecObject(const FieldEvaluateExecType& evaluator,
vtkm::FloatDefault stepLength,
vtkm::FloatDefault tolerance)
: Superclass(evaluator, stepLength, tolerance)
{
}
VTKM_EXEC
IntegratorStatus CheckStep(vtkm::Particle* particle,
vtkm::FloatDefault stepLength,
vtkm::Vec3f& velocity) const
{
auto time = particle->Time;
auto inpos = particle->Pos;
vtkm::VecVariable<vtkm::Vec3f, 2> vectors;
GridEvaluatorStatus status = this->Evaluator.Evaluate(inpos, time, vectors);
velocity = particle->Velocity(vectors, stepLength);
return IntegratorStatus(status);
}
};
private:
FieldEvaluateType Evaluator;
protected:
VTKM_CONT virtual void PrepareForExecutionImpl(
vtkm::cont::DeviceAdapterId device,
vtkm::cont::VirtualObjectHandle<IntegratorBase::ExecObject>& execObjectHandle,
vtkm::cont::Token& token) const override
{
vtkm::cont::TryExecuteOnDevice(device,
detail::IntegratorPrepareForExecutionFunctor<ExecObject>(),
execObjectHandle,
this->Evaluator,
this->StepLength,
this->Tolerance,
token);
}
}; //EulerIntegrator
} //namespace particleadvection
} //namespace worklet
} //namespace vtkm
#endif // vtk_m_worklet_particleadvection_EulerIntegrator_h

@ -24,7 +24,7 @@
#include <vtkm/cont/DeviceAdapter.h>
#include <vtkm/worklet/particleadvection/CellInterpolationHelper.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
namespace vtkm
{

@ -26,7 +26,7 @@
#include <vtkm/worklet/particleadvection/CellInterpolationHelper.h>
#include <vtkm/worklet/particleadvection/Field.h>
#include <vtkm/worklet/particleadvection/GridEvaluatorStatus.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
namespace vtkm
{

@ -10,8 +10,8 @@
//
//=============================================================================
#ifndef vtk_m_worklet_particleadvection_Integrators_h
#define vtk_m_worklet_particleadvection_Integrators_h
#ifndef vtk_m_worklet_particleadvection_IntegratorBase_h
#define vtk_m_worklet_particleadvection_IntegratorBase_h
#include <limits>
@ -34,14 +34,14 @@ namespace worklet
namespace particleadvection
{
class Integrator : public vtkm::cont::ExecutionObjectBase
class IntegratorBase : public vtkm::cont::ExecutionObjectBase
{
protected:
VTKM_CONT
Integrator() = default;
IntegratorBase() = default;
VTKM_CONT
Integrator(vtkm::FloatDefault stepLength)
IntegratorBase(vtkm::FloatDefault stepLength)
: StepLength(stepLength)
{
}
@ -244,7 +244,7 @@ struct IntegratorPrepareForExecutionFunctor
template <typename Device, typename EvaluatorType>
VTKM_CONT bool operator()(
Device,
vtkm::cont::VirtualObjectHandle<Integrator::ExecObject>& execObjectHandle,
vtkm::cont::VirtualObjectHandle<IntegratorBase::ExecObject>& execObjectHandle,
const EvaluatorType& evaluator,
vtkm::FloatDefault stepLength,
vtkm::FloatDefault tolerance,
@ -257,185 +257,9 @@ struct IntegratorPrepareForExecutionFunctor
}
};
} // namespace detail
template <typename FieldEvaluateType>
class RK4Integrator : public Integrator
{
public:
VTKM_CONT
RK4Integrator() = default;
VTKM_CONT
RK4Integrator(const FieldEvaluateType& evaluator, vtkm::FloatDefault stepLength)
: Integrator(stepLength)
, Evaluator(evaluator)
{
}
template <typename Device>
class ExecObject : public Integrator::ExecObjectBaseImpl<
vtkm::cont::internal::ExecutionObjectType<FieldEvaluateType, Device>,
typename RK4Integrator::template ExecObject<Device>>
{
VTKM_IS_DEVICE_ADAPTER_TAG(Device);
using FieldEvaluateExecType =
vtkm::cont::internal::ExecutionObjectType<FieldEvaluateType, Device>;
using Superclass =
Integrator::ExecObjectBaseImpl<FieldEvaluateExecType,
typename RK4Integrator::template ExecObject<Device>>;
public:
VTKM_EXEC_CONT
ExecObject(const FieldEvaluateExecType& evaluator,
vtkm::FloatDefault stepLength,
vtkm::FloatDefault tolerance)
: Superclass(evaluator, stepLength, tolerance)
{
}
VTKM_EXEC
IntegratorStatus CheckStep(vtkm::Particle* particle,
vtkm::FloatDefault stepLength,
vtkm::Vec3f& velocity) const
{
auto time = particle->Time;
auto inpos = particle->Pos;
vtkm::FloatDefault boundary = this->Evaluator.GetTemporalBoundary(static_cast<vtkm::Id>(1));
if ((time + stepLength + vtkm::Epsilon<vtkm::FloatDefault>() - 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<vtkm::FloatDefault>(2));
vtkm::FloatDefault var2 = time + var1;
vtkm::FloatDefault var3 = time + stepLength;
vtkm::Vec3f v1 = vtkm::TypeTraits<vtkm::Vec3f>::ZeroInitialization();
vtkm::Vec3f v2 = v1, v3 = v1, v4 = v1;
vtkm::VecVariable<vtkm::Vec3f, 2> k1, k2, k3, k4;
GridEvaluatorStatus evalStatus;
evalStatus = this->Evaluator.Evaluate(inpos, time, k1);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
v1 = particle->Velocity(k1, stepLength);
evalStatus = this->Evaluator.Evaluate(inpos + var1 * v1, var2, k2);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
v2 = particle->Velocity(k2, stepLength);
evalStatus = this->Evaluator.Evaluate(inpos + var1 * v2, var2, k3);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
v3 = particle->Velocity(k3, stepLength);
evalStatus = this->Evaluator.Evaluate(inpos + stepLength * v3, var3, k4);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
v4 = particle->Velocity(k4, stepLength);
velocity = (v1 + 2 * v2 + 2 * v3 + v4) / static_cast<vtkm::FloatDefault>(6);
return IntegratorStatus(true, false, evalStatus.CheckTemporalBounds());
}
};
private:
FieldEvaluateType Evaluator;
protected:
VTKM_CONT virtual void PrepareForExecutionImpl(
vtkm::cont::DeviceAdapterId device,
vtkm::cont::VirtualObjectHandle<Integrator::ExecObject>& execObjectHandle,
vtkm::cont::Token& token) const override
{
vtkm::cont::TryExecuteOnDevice(device,
detail::IntegratorPrepareForExecutionFunctor<ExecObject>(),
execObjectHandle,
this->Evaluator,
this->StepLength,
this->Tolerance,
token);
}
};
template <typename FieldEvaluateType>
class EulerIntegrator : public Integrator
{
public:
EulerIntegrator() = default;
VTKM_CONT
EulerIntegrator(const FieldEvaluateType& evaluator, const vtkm::FloatDefault stepLength)
: Integrator(stepLength)
, Evaluator(evaluator)
{
}
template <typename Device>
class ExecObject : public Integrator::ExecObjectBaseImpl<
vtkm::cont::internal::ExecutionObjectType<FieldEvaluateType, Device>,
typename EulerIntegrator::template ExecObject<Device>>
{
VTKM_IS_DEVICE_ADAPTER_TAG(Device);
using FieldEvaluateExecType =
vtkm::cont::internal::ExecutionObjectType<FieldEvaluateType, Device>;
using Superclass =
Integrator::ExecObjectBaseImpl<FieldEvaluateExecType,
typename EulerIntegrator::template ExecObject<Device>>;
public:
VTKM_EXEC_CONT
ExecObject(const FieldEvaluateExecType& evaluator,
vtkm::FloatDefault stepLength,
vtkm::FloatDefault tolerance)
: Superclass(evaluator, stepLength, tolerance)
{
}
VTKM_EXEC
IntegratorStatus CheckStep(vtkm::Particle* particle,
vtkm::FloatDefault stepLength,
vtkm::Vec3f& velocity) const
{
auto time = particle->Time;
auto inpos = particle->Pos;
vtkm::VecVariable<vtkm::Vec3f, 2> vectors;
GridEvaluatorStatus status = this->Evaluator.Evaluate(inpos, time, vectors);
velocity = particle->Velocity(vectors, stepLength);
return IntegratorStatus(status);
}
};
private:
FieldEvaluateType Evaluator;
protected:
VTKM_CONT virtual void PrepareForExecutionImpl(
vtkm::cont::DeviceAdapterId device,
vtkm::cont::VirtualObjectHandle<Integrator::ExecObject>& execObjectHandle,
vtkm::cont::Token& token) const override
{
vtkm::cont::TryExecuteOnDevice(device,
detail::IntegratorPrepareForExecutionFunctor<ExecObject>(),
execObjectHandle,
this->Evaluator,
this->StepLength,
this->Tolerance,
token);
}
}; //EulerIntegrator
} //namespace detail
} //namespace particleadvection
} //namespace worklet
} //namespace vtkm
#endif // vtk_m_worklet_particleadvection_Integrators_h
#endif // vtk_m_worklet_particleadvection_IntegratorBase_h

@ -22,9 +22,8 @@
#include <vtkm/Particle.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
#include <vtkm/worklet/particleadvection/Particles.h>
//#include <vtkm/worklet/particleadvection/ParticlesAOS.h>
#ifdef VTKM_CUDA
#include <vtkm/cont/cuda/internal/ScopedCudaStackSize.h>

@ -0,0 +1,137 @@
//=============================================================================
//
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
//=============================================================================
#ifndef vtk_m_worklet_particleadvection_RK4Integrator_h
#define vtk_m_worklet_particleadvection_RK4Integrator_h
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
namespace vtkm
{
namespace worklet
{
namespace particleadvection
{
template <typename FieldEvaluateType>
class RK4Integrator : public IntegratorBase
{
public:
VTKM_CONT
RK4Integrator() = default;
VTKM_CONT
RK4Integrator(const FieldEvaluateType& evaluator, vtkm::FloatDefault stepLength)
: IntegratorBase(stepLength)
, Evaluator(evaluator)
{
}
template <typename Device>
class ExecObject : public IntegratorBase::ExecObjectBaseImpl<
vtkm::cont::internal::ExecutionObjectType<FieldEvaluateType, Device>,
typename RK4Integrator::template ExecObject<Device>>
{
VTKM_IS_DEVICE_ADAPTER_TAG(Device);
using FieldEvaluateExecType =
vtkm::cont::internal::ExecutionObjectType<FieldEvaluateType, Device>;
using Superclass =
IntegratorBase::ExecObjectBaseImpl<FieldEvaluateExecType,
typename RK4Integrator::template ExecObject<Device>>;
public:
VTKM_EXEC_CONT
ExecObject(const FieldEvaluateExecType& evaluator,
vtkm::FloatDefault stepLength,
vtkm::FloatDefault tolerance)
: Superclass(evaluator, stepLength, tolerance)
{
}
VTKM_EXEC
IntegratorStatus CheckStep(vtkm::Particle* particle,
vtkm::FloatDefault stepLength,
vtkm::Vec3f& velocity) const
{
auto time = particle->Time;
auto inpos = particle->Pos;
vtkm::FloatDefault boundary = this->Evaluator.GetTemporalBoundary(static_cast<vtkm::Id>(1));
if ((time + stepLength + vtkm::Epsilon<vtkm::FloatDefault>() - 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<vtkm::FloatDefault>(2));
vtkm::FloatDefault var2 = time + var1;
vtkm::FloatDefault var3 = time + stepLength;
vtkm::Vec3f v1 = vtkm::TypeTraits<vtkm::Vec3f>::ZeroInitialization();
vtkm::Vec3f v2 = v1, v3 = v1, v4 = v1;
vtkm::VecVariable<vtkm::Vec3f, 2> k1, k2, k3, k4;
GridEvaluatorStatus evalStatus;
evalStatus = this->Evaluator.Evaluate(inpos, time, k1);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
v1 = particle->Velocity(k1, stepLength);
evalStatus = this->Evaluator.Evaluate(inpos + var1 * v1, var2, k2);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
v2 = particle->Velocity(k2, stepLength);
evalStatus = this->Evaluator.Evaluate(inpos + var1 * v2, var2, k3);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
v3 = particle->Velocity(k3, stepLength);
evalStatus = this->Evaluator.Evaluate(inpos + stepLength * v3, var3, k4);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
v4 = particle->Velocity(k4, stepLength);
velocity = (v1 + 2 * v2 + 2 * v3 + v4) / static_cast<vtkm::FloatDefault>(6);
return IntegratorStatus(true, false, evalStatus.CheckTemporalBounds());
}
};
private:
FieldEvaluateType Evaluator;
protected:
VTKM_CONT virtual void PrepareForExecutionImpl(
vtkm::cont::DeviceAdapterId device,
vtkm::cont::VirtualObjectHandle<IntegratorBase::ExecObject>& execObjectHandle,
vtkm::cont::Token& token) const override
{
vtkm::cont::TryExecuteOnDevice(device,
detail::IntegratorPrepareForExecutionFunctor<ExecObject>(),
execObjectHandle,
this->Evaluator,
this->StepLength,
this->Tolerance,
token);
}
};
} //namespace particleadvection
} //namespace worklet
} //namespace vtkm
#endif // vtk_m_worklet_particleadvection_RK4Integrator_h

@ -18,10 +18,12 @@
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/io/VTKDataSetReader.h>
#include <vtkm/worklet/ParticleAdvection.h>
#include <vtkm/worklet/particleadvection/EulerIntegrator.h>
#include <vtkm/worklet/particleadvection/Field.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
#include <vtkm/worklet/particleadvection/Particles.h>
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
#include <random>

@ -19,7 +19,7 @@
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/worklet/ParticleAdvection.h>
#include <vtkm/worklet/particleadvection/Field.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
#include <vtkm/worklet/particleadvection/Particles.h>
#include <vtkm/worklet/particleadvection/TemporalGridEvaluators.h>