Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-10-05 01:49:02 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Removing XGC Filter

Browse Source
This commit is contained in:
Abhishek Yenpure 2023-07-11 13:20:05 -07:00 committed by Abhishek Yenpure
parent 9d840a88bd
commit 9c83bd983b
9 changed files with 0 additions and 1393 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
vtkm/filter/flow/CMakeLists.txt
View File

@ -21,7 +21,6 @@ set(flow_headers
Streamline.h
StreamSurface.h
WarpXStreamline.h
XGCPoincare.h
)
set(flow_sources
@ -42,7 +41,6 @@ set(flow_device_sources
PathParticle.cxx
Streamline.cxx
WarpXStreamline.cxx
XGCPoincare.cxx
)
vtkm_library(

90
vtkm/filter/flow/XGCPoincare.cxx
View File

@ -1,90 +0,0 @@
//============================================================================
// 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.
//============================================================================
#include <vtkm/filter/flow/XGCPoincare.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
VTKM_CONT XGCPoincare::FieldType XGCPoincare::GetField(const vtkm::cont::DataSet& dataset) const
{
const auto& field1 = this->GetActiveFieldName(0);
const auto& field2 = this->GetActiveFieldName(1);
const auto& field3 = this->GetActiveFieldName(2);
const auto& field4 = this->GetActiveFieldName(3);
const auto& field5 = this->GetActiveFieldName(4);
const auto& field6 = this->GetActiveFieldName(5);
/*
XGCField(const FieldComponentType& as_ff,
const FieldVecType& _dAs_ff_rzp,
const FieldComponentType& coeff_1D,
const FieldComponentType& coeff_2D,
const FieldComponentType& psi,
const FieldVecType& _B_rzp,
const XGCParams& params,
bool useDeltaBScale,
bool useBScale,
const Association assoc)
*/
CompArrayType as_ff;
CompArrayType coeff_1D;
CompArrayType coeff_2D;
CompArrayType psi;
VecArrayType dAs_ff_rzp;
VecArrayType b_rzp;
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(field1).GetData(), as_ff);
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(field2).GetData(), dAs_ff_rzp);
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(field3).GetData(), coeff_1D);
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(field4).GetData(), coeff_2D);
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(field5).GetData(), b_rzp);
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(field6).GetData(), psi);
vtkm::worklet::flow::XGCParams params = this->GetXGCParams();
vtkm::FloatDefault period = this->GetPeriod();
return XGCPoincare::FieldType(
as_ff, dAs_ff_rzp, coeff_1D, coeff_2D, b_rzp, psi, params, period, false, 1.0, false, 1.0);
}
VTKM_CONT XGCPoincare::TerminationType XGCPoincare::GetTermination(
const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return XGCPoincare::TerminationType(this->NumberOfSteps, this->MaxPunctures);
}
VTKM_CONT XGCPoincare::AnalysisType XGCPoincare::GetAnalysis(
const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
CompArrayType coeff_1D;
CompArrayType coeff_2D;
const auto& field3 = this->GetActiveFieldName(2);
const auto& field4 = this->GetActiveFieldName(3);
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(field3).GetData(), coeff_1D);
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(field4).GetData(), coeff_2D);
vtkm::worklet::flow::XGCParams params = this->GetXGCParams();
vtkm::FloatDefault period = this->GetPeriod();
return XGCPoincare::AnalysisType(
this->NumberOfSteps, this->MaxPunctures, coeff_1D, coeff_2D, params, period);
}
}
}
} // namespace vtkm::filter::flow

106
vtkm/filter/flow/XGCPoincare.h
View File

@ -1,106 +0,0 @@
//============================================================================
// 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_filter_flow_XGCPoincare_h
#define vtk_m_filter_flow_XGCPoincare_h
#include <vtkm/filter/flow/FilterParticleAdvectionSteadyState.h>
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
#include <vtkm/filter/flow/worklet/Analysis.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/Termination.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
class XGCPoincare;
template <>
struct FlowTraits<XGCPoincare>
{
using ParticleType = vtkm::Particle;
using TerminationType = vtkm::worklet::flow::PoincareTermination;
using AnalysisType = vtkm::worklet::flow::XGCPoincare<ParticleType>;
using CompArrayType = vtkm::cont::ArrayHandle<vtkm::FloatDefault>;
using VecArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::flow::XGCField<CompArrayType, VecArrayType>;
};
/// \brief Advect particles in a vector field.
/// Takes as input a vector field and seed locations and generates the
/// end points for each seed through the vector field.
class VTKM_FILTER_FLOW_EXPORT XGCPoincare
: public vtkm::filter::flow::FilterParticleAdvectionSteadyState<XGCPoincare>
{
public:
using ParticleType = typename FlowTraits<XGCPoincare>::ParticleType;
using TerminationType = typename FlowTraits<XGCPoincare>::TerminationType;
using AnalysisType = typename FlowTraits<XGCPoincare>::AnalysisType;
using CompArrayType = typename FlowTraits<XGCPoincare>::CompArrayType;
using VecArrayType = typename FlowTraits<XGCPoincare>::VecArrayType;
using FieldType = typename FlowTraits<XGCPoincare>::FieldType;
VTKM_CONT FieldType GetField(const vtkm::cont::DataSet& data) const;
VTKM_CONT TerminationType GetTermination(const vtkm::cont::DataSet& data) const;
VTKM_CONT AnalysisType GetAnalysis(const vtkm::cont::DataSet& data) const;
VTKM_CONT void SetXGCParams(const vtkm::worklet::flow::XGCParams& params)
{
this->Params = params;
}
VTKM_CONT vtkm::worklet::flow::XGCParams GetXGCParams() const { return this->Params; }
VTKM_CONT void SetPeriod(vtkm::FloatDefault period) { this->Period = period; }
VTKM_CONT vtkm::FloatDefault GetPeriod() const { return this->Period; }
VTKM_CONT void SetMaxPunctures(vtkm::Id maxPunctures) { this->MaxPunctures = maxPunctures; }
VTKM_CONT vtkm::Id GetMaxPunctures() const { return this->MaxPunctures; }
// As_ff
VTKM_CONT void SetField1(const std::string& name) { this->SetActiveField(0, name); }
// dAs_ff
VTKM_CONT void SetField2(const std::string& name) { this->SetActiveField(1, name); }
// Bfield
VTKM_CONT void SetField3(const std::string& name) { this->SetActiveField(2, name); }
// Psi
VTKM_CONT void SetField4(const std::string& name) { this->SetActiveField(3, name); }
// coeff_1d
VTKM_CONT void SetField5(const std::string& name) { this->SetActiveField(4, name); }
// coeff_2d
VTKM_CONT void SetField6(const std::string& name) { this->SetActiveField(5, name); }
VTKM_CONT std::string GetField1() const { return this->GetActiveFieldName(0); }
VTKM_CONT std::string GetField2() const { return this->GetActiveFieldName(1); }
VTKM_CONT std::string GetField3() const { return this->GetActiveFieldName(2); }
VTKM_CONT std::string GetField4() const { return this->GetActiveFieldName(3); }
VTKM_CONT std::string GetField5() const { return this->GetActiveFieldName(4); }
VTKM_CONT std::string GetField6() const { return this->GetActiveFieldName(5); }
private:
vtkm::worklet::flow::XGCParams Params;
vtkm::FloatDefault Period;
vtkm::Id MaxPunctures;
};
}
}
} // namespace vtkm::filter::flow
#endif // vtk_m_filter_flow_XGCPoincare_h

2
vtkm/filter/flow/worklet/Analysis.h
View File

@ -407,6 +407,4 @@ private:
} // namespace worklet
} // namespace vtkm
#include <vtkm/filter/flow/worklet/XGCPoincare.h>
#endif //vtkm_worklet_particleadvection_analysis

3
vtkm/filter/flow/worklet/CMakeLists.txt
View File

@ -26,9 +26,6 @@ set(headers
StreamSurface.h
TemporalGridEvaluators.h
Termination.h
XGCHelper.h
XGCField.h
XGCPoincare.h
)
#-----------------------------------------------------------------------------

2
vtkm/filter/flow/worklet/Field.h
View File

@ -237,6 +237,4 @@ private:
}
} //vtkm::worklet::flow
#include <vtkm/filter/flow/worklet/XGCField.h>
#endif //vtkm_filter_flow_worklet_Field_h

400
vtkm/filter/flow/worklet/XGCField.h
View File

@ -1,400 +0,0 @@
//============================================================================
// 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 vtkm_filter_flow_worklet_XGCField_h
#define vtkm_filter_flow_worklet_XGCField_h
#include <vtkm/Geometry.h>
#include <vtkm/Matrix.h>
#include <vtkm/Types.h>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/VecVariable.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ExecutionObjectBase.h>
#include <vtkm/exec/CellInterpolate.h>
#include <vtkm/filter/flow/worklet/XGCHelper.h>
namespace vtkm
{
namespace worklet
{
namespace flow
{
template <typename FieldComponentType, typename FieldVecType>
class ExecutionXGCField
{
public:
using ComponentPortalType = typename FieldComponentType::ReadPortalType;
using VecPortalType = typename FieldVecType::ReadPortalType;
using Association = vtkm::cont::Field::Association;
using DelegateToField = std::true_type;
using Ray3f = vtkm::Ray<vtkm::FloatDefault, 3, true>;
VTKM_CONT
ExecutionXGCField(const FieldComponentType& as_ff,
const FieldVecType& _dAs_ff_rzp,
const FieldComponentType& coeff_1D,
const FieldComponentType& coeff_2D,
const FieldVecType& _B_rzp,
const FieldComponentType& psi,
const XGCParams& params,
vtkm::FloatDefault period,
bool useDeltaBScale,
vtkm::FloatDefault deltaBScale,
bool useBScale,
vtkm::FloatDefault bScale,
const Association assoc,
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token)
: As_ff(as_ff.PrepareForInput(device, token))
, dAs_ff_rzp(_dAs_ff_rzp.PrepareForInput(device, token))
, Coeff_1D(coeff_1D.PrepareForInput(device, token))
, Coeff_2D(coeff_2D.PrepareForInput(device, token))
, Psi(psi.PrepareForInput(device, token))
, B_rzp(_B_rzp.PrepareForInput(device, token))
, Params(params)
, Period(period)
, UseDeltaBScale(useDeltaBScale)
, DeltaBScale(deltaBScale)
, UseBScale(useBScale)
, BScale(bScale)
, Assoc(assoc)
{
}
private:
template <typename LocatorType, typename InterpolationHelper>
VTKM_EXEC bool PtLoc(const vtkm::Vec3f& ptRZ,
const LocatorType& locator,
const InterpolationHelper& helper,
vtkm::Vec3f& param,
vtkm::Vec<vtkm::Id, 3>& vIds) const
{
vtkm::Id cellId;
vtkm::ErrorCode status;
status = locator.FindCell(ptRZ, cellId, param);
if (status != vtkm::ErrorCode::Success)
{
return false;
}
vtkm::UInt8 cellShape;
vtkm::IdComponent nVerts;
vtkm::VecVariable<vtkm::Id, 8> ptIndices;
helper.GetCellInfo(cellId, cellShape, nVerts, ptIndices);
vIds[0] = ptIndices[0];
vIds[1] = ptIndices[1];
vIds[2] = ptIndices[2];
return true;
}
public:
VTKM_EXEC Association GetAssociation() const { return this->Assoc; }
VTKM_EXEC void GetValue(const vtkm::Id vtkmNotUsed(cellId),
vtkm::VecVariable<vtkm::Vec3f, 2>& vtkmNotUsed(value)) const
{
//TODO Raise Error : XGC Field should not allow this path
}
VTKM_EXEC void GetValue(const vtkm::VecVariable<vtkm::Id, 8>& vtkmNotUsed(indices),
const vtkm::Id vtkmNotUsed(vertices),
const vtkm::Vec3f& vtkmNotUsed(parametric),
const vtkm::UInt8 vtkmNotUsed(cellShape),
vtkm::VecVariable<vtkm::Vec3f, 2>& vtkmNotUsed(value)) const
{
//TODO Raise Error : XGC Field should not allow this path
}
template <typename Point, typename Locator, typename Helper>
VTKM_EXEC bool GetValue(const Point& point,
const vtkm::FloatDefault& time,
vtkm::VecVariable<Point, 2>& out,
const Locator& locator,
const Helper& helper) const
{
// not used
(void)time;
auto R = point[0];
auto Phi = point[1];
auto Z = point[2];
ParticleMetaData metadata;
if (!HighOrderB(point, metadata, this->Params, this->Coeff_1D, this->Coeff_2D))
return false;
vtkm::Id planeIdx0, planeIdx1, numRevs;
vtkm::FloatDefault phiN, Phi0, Phi1, T;
GetPlaneIdx(
Phi, this->Period, this->Params, phiN, planeIdx0, planeIdx1, Phi0, Phi1, numRevs, T);
vtkm::Vec3f B0_rpz(metadata.B0_rzp[0], metadata.B0_rzp[2], metadata.B0_rzp[1]);
vtkm::Vec3f ff_pt_rpz;
CalcFieldFollowingPt(
{ R, phiN, Z }, B0_rpz, Phi0, Phi1, this->Params, this->Coeff_1D, this->Coeff_2D, ff_pt_rpz);
//Now, interpolate between Phi_i and Phi_i+1
vtkm::FloatDefault T01 = (phiN - Phi0) / (Phi1 - Phi0);
vtkm::FloatDefault T10 = 1.0f - T01;
//Get vec at Phi0 and Phi1.
//x_ff is in rzp
//vtkm::Vec3f x_ff_rzp(ptOnMidPlane_rpz[0], ptOnMidPlane_rpz[2], 0);
vtkm::Vec3f x_ff_rzp(ff_pt_rpz[0], ff_pt_rpz[2], 0);
// Offsets into As and DAs to jump to the right plane.
int offsets[2];
offsets[0] = planeIdx0 * this->Params.NumNodes * 2;
offsets[1] = planeIdx0 * this->Params.NumNodes * 2 + this->Params.NumNodes;
const vtkm::FloatDefault basis = 0.0f;
//auto B0_R = B0_rzp[0];
//auto B0_Z = B0_rzp[1];
//auto x_ff_R = x_ff_rzp[0];
//auto x_ff_Z = x_ff_rzp[1];
//gradPsi: pt on mid plane? (question)
//dPsi/dR = B0_Z * R
//dPsi/dZ = -B0_R * R;
//vtkm::Vec3f gradPsi_rzp(B0_Z * x_ff_R, -B0_R * x_ff_R, 0);
//use high order...
vtkm::Vec3f gradPsi_rzp = metadata.gradPsi_rzp;
vtkm::FloatDefault gammaPsi = 1.0f / vtkm::Magnitude(gradPsi_rzp);
vtkm::Vec2f rvec(0, 0), zvec(0, 0);
rvec[0] = basis + (1.0 - basis) * gammaPsi * gradPsi_rzp[0];
rvec[1] = (1.0 - basis) * gammaPsi * gradPsi_rzp[1];
zvec[0] = (1.0 - basis) * gammaPsi * (-gradPsi_rzp[1]);
zvec[1] = basis + (1.0 - basis) * gammaPsi * gradPsi_rzp[0];
//Get the vectors in the ff coordinates.
//auto dAs_ff_rzp = EvalVector(ds, locator, {x_ff_rzp, x_ff_rzp}, "dAs_ff_rzp", offsets);
//auto dAs_ff0_rzp = dAs_ff_rzp[0];
//auto dAs_ff1_rzp = dAs_ff_rzp[1];
vtkm::Vec3f x_ff_param;
vtkm::Vec<vtkm::Id, 3> x_ff_vids;
//Hotspot..
if (!this->PtLoc(x_ff_rzp, locator, helper, x_ff_param, x_ff_vids))
return false;
// Eval actual values for DAsPhiFF from the triangular planes
auto dAs_ff0_rzp = Eval(this->dAs_ff_rzp, offsets[0], x_ff_param, x_ff_vids);
auto dAs_ff1_rzp = Eval(this->dAs_ff_rzp, offsets[1], x_ff_param, x_ff_vids);
vtkm::FloatDefault wphi[2] = { T10, T01 }; //{T01, T10};
vtkm::Vec3f gradAs_rpz;
//vec.r = wphi[0]*( rvec[0]*V.r[0] + zvec[0]*V.z[0]) +
// wphi[1]*( rvec[0]*V.r[1] + zvec[0]*v.z[1]);
//vec.p = wphi[0]*V.phi[0] +
// whpi[1]*V.phi[1];
//vec.z = wphi[0]*( rvec[1]*V.r[0] + zvec[1]*V.z[0]) +
// wphi[1]*( rvec[1]*V.r[1] + zvec[1]*V.Z[1]);
gradAs_rpz[0] = wphi[0] * (rvec[0] * dAs_ff0_rzp[0] + zvec[0] * dAs_ff0_rzp[1]) +
wphi[1] * (rvec[0] * dAs_ff1_rzp[0] + zvec[0] * dAs_ff1_rzp[1]);
gradAs_rpz[1] = wphi[0] * dAs_ff0_rzp[2] + wphi[1] * dAs_ff1_rzp[2];
gradAs_rpz[2] = wphi[0] * (rvec[1] * dAs_ff0_rzp[0] + zvec[1] * dAs_ff0_rzp[1]) +
wphi[1] * (rvec[1] * dAs_ff1_rzp[0] + zvec[1] * dAs_ff1_rzp[1]);
vtkm::FloatDefault BMag = vtkm::Magnitude(metadata.B0_rzp);
//project using bfield.
//gradAs.Phi = (gradAs.Phi * BMag - gradAs.R*B0_pos.R - gradAs.Z*B0_pos.Z) / B0_pos.Phi
gradAs_rpz[1] = (gradAs_rpz[1] * BMag - gradAs_rpz[0] * metadata.B0_rzp[0] -
gradAs_rpz[2] * metadata.B0_rzp[1]) /
metadata.B0_rzp[2];
//deltaB = AsCurl(bhat) + gradAs x bhat.
//std::vector<int> off = {planeIdx0*this->NumNodes};
//vtkm::Vec3f AsCurl_bhat_rzp = EvalVector(ds, locator, {x_ff_rzp}, "AsCurlBHat_RZP", off)[0];
//auto AsCurl_bhat_rzp = this->EvalV(AsCurlBHat_RZP, 0, x_ff_vids, x_ff_param);
//
auto As_ff0 = Eval(this->As_ff, offsets[0], x_ff_param, x_ff_vids);
auto As_ff1 = Eval(this->As_ff, offsets[1], x_ff_param, x_ff_vids);
// Interpolated value of As on plane 10 & 1.
// wPhi is distance of the plane from the point.
vtkm::FloatDefault As = wphi[0] * As_ff0 + wphi[1] * As_ff1;
auto AsCurl_bhat_rzp = As * metadata.curl_nb_rzp;
auto bhat_rzp = vtkm::Normal(metadata.B0_rzp);
vtkm::Vec3f gradAs_rzp(gradAs_rpz[0], gradAs_rpz[2], gradAs_rpz[1]);
vtkm::Vec3f deltaB_rzp = AsCurl_bhat_rzp + vtkm::Cross(gradAs_rzp, bhat_rzp);
if (this->UseDeltaBScale)
deltaB_rzp = deltaB_rzp * this->DeltaBScale;
vtkm::Vec3f B0_rzp = metadata.B0_rzp;
if (this->UseBScale)
B0_rzp = B0_rzp * this->BScale;
deltaB_rzp[2] /= R;
B0_rzp[2] /= R;
vtkm::Vec3f vec_rzp = B0_rzp + deltaB_rzp;
vtkm::Vec3f vec_rpz(vec_rzp[0], vec_rzp[2], vec_rzp[1]);
out = vtkm::make_Vec(vec_rpz);
return true;
}
private:
ComponentPortalType As_ff;
VecPortalType dAs_ff_rzp;
ComponentPortalType Coeff_1D;
ComponentPortalType Coeff_2D;
ComponentPortalType Psi;
VecPortalType B_rzp;
vtkm::FloatDefault Period;
bool UseDeltaBScale;
vtkm::FloatDefault DeltaBScale = 1.0;
bool UseBScale;
vtkm::FloatDefault BScale = 1.0;
Association Assoc;
XGCParams Params;
};
template <typename FieldComponentType, typename FieldVecType>
class XGCField : public vtkm::cont::ExecutionObjectBase
{
public:
using ExecutionType = ExecutionXGCField<FieldComponentType, FieldVecType>;
using Association = vtkm::cont::Field::Association;
VTKM_CONT
XGCField() = default;
VTKM_CONT
XGCField(const FieldComponentType& as_ff,
const FieldVecType& _dAs_ff_rzp,
const FieldComponentType& coeff_1D,
const FieldComponentType& coeff_2D,
const FieldVecType& _B_rzp,
const FieldComponentType& psi,
const XGCParams& params,
vtkm::FloatDefault period,
bool useDeltaBScale,
vtkm::FloatDefault deltaBScale,
bool useBScale,
vtkm::FloatDefault bScale)
: As_ff(as_ff)
, dAs_ff_rzp(_dAs_ff_rzp)
, Coeff_1D(coeff_1D)
, Coeff_2D(coeff_2D)
, Psi(psi)
, B_rzp(_B_rzp)
, Params(params)
, Period(period)
, UseDeltaBScale(useDeltaBScale)
, DeltaBScale(deltaBScale)
, UseBScale(useBScale)
, BScale(bScale)
, Assoc(vtkm::cont::Field::Association::Points)
{
}
VTKM_CONT
XGCField(const FieldComponentType& as_ff,
const FieldVecType& _dAs_ff_rzp,
const FieldComponentType& coeff_1D,
const FieldComponentType& coeff_2D,
const FieldVecType& _B_rzp,
const FieldComponentType& psi,
const XGCParams& params,
vtkm::FloatDefault period,
bool useDeltaBScale,
vtkm::FloatDefault deltaBScale,
bool useBScale,
vtkm::FloatDefault bScale,
const Association assoc)
: As_ff(as_ff)
, dAs_ff_rzp(_dAs_ff_rzp)
, Coeff_1D(coeff_1D)
, Coeff_2D(coeff_2D)
, Psi(psi)
, B_rzp(_B_rzp)
, Params(params)
, Period(period)
, UseDeltaBScale(useDeltaBScale)
, DeltaBScale(deltaBScale)
, UseBScale(useBScale)
, BScale(bScale)
, Assoc(assoc)
{
if (assoc != Association::Points && assoc != Association::Cells)
throw("Unsupported field association");
}
VTKM_CONT
const ExecutionType PrepareForExecution(vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
{
return ExecutionType(this->As_ff,
this->dAs_ff_rzp,
this->Coeff_1D,
this->Coeff_2D,
this->B_rzp,
this->Psi,
this->Params,
this->Period,
this->UseDeltaBScale,
this->DeltaBScale,
this->UseBScale,
this->BScale,
this->Assoc,
device,
token);
}
private:
// For turbulance in the magnetic field (B is constant, but this emulates deltaB)
FieldComponentType As_ff;
FieldVecType dAs_ff_rzp;
// For Higher order interpolation
// on RZ uniform grid, not the triangle grid.
FieldComponentType Coeff_1D;
FieldComponentType Coeff_2D;
// Constant magnetic field on the triangle mesh
FieldComponentType Psi;
// B field at each vertex of the triangular field
FieldVecType B_rzp;
XGCParams Params;
vtkm::FloatDefault Period;
bool UseDeltaBScale;
vtkm::FloatDefault DeltaBScale = 1.0;
bool UseBScale;
vtkm::FloatDefault BScale = 1.0;
Association Assoc;
};
}
}
} //vtkm::worklet::flow
#endif //vtkm_filter_flow_worklet_XGCField_h

507
vtkm/filter/flow/worklet/XGCHelper.h
View File

@ -1,507 +0,0 @@
//============================================================================
// 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 vtkm_filter_flow_worklet_XGCHelper_h
#define vtkm_filter_flow_worklet_XGCHelper_h
#include <vtkm/Geometry.h>
#include <vtkm/Matrix.h>
namespace vtkm
{
namespace worklet
{
namespace flow
{
struct XGCParams
{
vtkm::Id NumPlanes = -1;
vtkm::Id NumNodes = -1;
vtkm::Id NumTri = -1;
vtkm::FloatDefault Period;
vtkm::FloatDefault dPhi;
int ncoeff_r, ncoeff_z, ncoeff_psi;
vtkm::FloatDefault min_psi, max_psi;
vtkm::FloatDefault sml_bp_sign = -1.0f;
vtkm::FloatDefault eq_axis_r = -1, eq_axis_z = -1, eq_x_psi = -1, eq_x_r = -1, eq_x_z = -1;
vtkm::FloatDefault eq_min_r = -1, eq_max_r = -1;
vtkm::FloatDefault eq_min_z = -1, eq_max_z = 1;
vtkm::FloatDefault psi_min = -1, psi_max = -1;
vtkm::Id eq_mr = -1, eq_mz = -1;
vtkm::Id eq_mpsi = -1;
vtkm::Id sml_wedge_n = -1;
vtkm::FloatDefault itp_min_psi = -1, itp_max_psi = -1;
vtkm::FloatDefault one_d_cub_dpsi_inv;
int nr, nz;
vtkm::FloatDefault rmin, zmin, rmax, zmax;
vtkm::FloatDefault dr, dz, dr_inv, dz_inv;
vtkm::FloatDefault eq_x2_z, eq_x2_r, eq_x2_psi;
vtkm::FloatDefault eq_x_slope, eq_x2_slope;
};
struct ParticleMetaData
{
ParticleMetaData() = default;
ParticleMetaData(const ParticleMetaData&) = default;
ParticleMetaData& operator=(const ParticleMetaData&) = default;
vtkm::Id3 PrevCell = { -1, -1, -1 };
vtkm::FloatDefault Psi = 0;
vtkm::FloatDefault dpsi_dr = 0, dpsi_dz = 0, d2psi_drdz = 0, d2psi_d2r = 0, d2psi_d2z = 0;
vtkm::Vec3f gradPsi_rzp, B0_rzp, curlB_rzp, curl_nb_rzp;
};
VTKM_EXEC
inline vtkm::Id GetIndex(const vtkm::FloatDefault& x,
const vtkm::Id& nx,
const vtkm::FloatDefault& xmin,
const vtkm::FloatDefault& dx_inv)
{
vtkm::Id idx = vtkm::Max(vtkm::Id(0), vtkm::Min(nx - 1, vtkm::Id((x - xmin) * dx_inv)));
return idx;
}
template <typename ScalarPortal>
VTKM_EXEC void EvalBicub2(const vtkm::FloatDefault& x,
const vtkm::FloatDefault& y,
const vtkm::FloatDefault& xc,
const vtkm::FloatDefault& yc,
const vtkm::Id& offset,
const ScalarPortal& Coeff_2D,
vtkm::FloatDefault& f00,
vtkm::FloatDefault& f10,
vtkm::FloatDefault& f01,
vtkm::FloatDefault& f11,
vtkm::FloatDefault& f20,
vtkm::FloatDefault& f02)
{
vtkm::FloatDefault dx = x - xc;
vtkm::FloatDefault dy = y - yc;
//fortran code.
f00 = f01 = f10 = f11 = f20 = f02 = 0.0f;
vtkm::FloatDefault xv[4] = { 1, dx, dx * dx, dx * dx * dx };
vtkm::FloatDefault yv[4] = { 1, dy, dy * dy, dy * dy * dy };
vtkm::FloatDefault fx[4] = { 0, 0, 0, 0 };
vtkm::FloatDefault dfx[4] = { 0, 0, 0, 0 };
vtkm::FloatDefault dfy[4] = { 0, 0, 0, 0 };
vtkm::FloatDefault dfx2[4] = { 0, 0, 0, 0 };
vtkm::FloatDefault dfy2[4] = { 0, 0, 0, 0 };
for (int j = 0; j < 4; j++)
{
for (int i = 0; i < 4; i++)
fx[j] = fx[j] + xv[i] * Coeff_2D.Get(offset + j * 4 + i); //acoeff[i][j];
for (int i = 1; i < 4; i++)
dfx[j] = dfx[j] +
vtkm::FloatDefault(i) * xv[i - 1] * Coeff_2D.Get(offset + j * 4 + i); //acoeff[i][j];
for (int i = 2; i < 4; i++)
dfx2[j] = dfx2[j] +
vtkm::FloatDefault(i * (i - 1)) * xv[i - 2] *
Coeff_2D.Get(offset + j * 4 + i); //acoeff[i][j];
}
for (int j = 0; j < 4; j++)
{
f00 = f00 + fx[j] * yv[j];
f10 = f10 + dfx[j] * yv[j];
f20 = f20 + dfx2[j] * yv[j];
}
for (int j = 1; j < 4; j++)
{
dfy[j] = vtkm::FloatDefault(j) * yv[j - 1];
f01 = f01 + fx[j] * dfy[j];
f11 = f11 + dfx[j] * dfy[j];
}
for (int j = 2; j < 4; j++)
{
dfy2[j] = vtkm::FloatDefault(j * (j - 1)) * yv[j - 2];
f02 = f02 + fx[j] * dfy2[j];
}
}
template <typename ScalarPortal>
VTKM_EXEC vtkm::FloatDefault I_interpol(const vtkm::FloatDefault& in_psi,
const vtkm::Id& ideriv,
const vtkm::Id& region,
const XGCParams& Params,
const ScalarPortal& Coeff_1D)
{
vtkm::FloatDefault psi;
if (region == 3)
{
psi = vtkm::Min(Params.eq_x_psi, Params.itp_max_psi);
if (ideriv != 0)
return 0.0;
}
else
{
psi = vtkm::Min(in_psi, Params.itp_max_psi);
psi = vtkm::Max(psi, Params.itp_min_psi);
}
vtkm::FloatDefault pn = psi * Params.one_d_cub_dpsi_inv;
vtkm::Id ip = floor(pn);
ip = vtkm::Min(vtkm::Max(ip, vtkm::Id(0)), vtkm::Id(Params.ncoeff_psi - 1));
vtkm::FloatDefault wp = pn - (vtkm::FloatDefault)(ip);
int idx = ip * 4;
//vtkm::FloatDefault acoef[4];
//acoef[0] = one_d_cub_acoef(ip).coeff[0];
//acoef[1] = one_d_cub_acoef(ip).coeff[1];
//acoef[2] = one_d_cub_acoef(ip).coeff[2];
//acoef[3] = one_d_cub_acoef(ip).coeff[3];
const vtkm::FloatDefault acoef[4] = {
Coeff_1D.Get(idx + 0), Coeff_1D.Get(idx + 1), Coeff_1D.Get(idx + 2), Coeff_1D.Get(idx + 3)
};
vtkm::FloatDefault iVal = 0.0;
if (ideriv == 0)
iVal = acoef[0] + (acoef[1] + (acoef[2] + acoef[3] * wp) * wp) * wp;
else if (ideriv == 1)
iVal = (acoef[1] + (2.0 * acoef[2] + 3.0 * acoef[3] * wp) * wp) * Params.one_d_cub_dpsi_inv;
return iVal * Params.sml_bp_sign;
}
VTKM_EXEC
inline bool IsRegion12(const vtkm::FloatDefault& R,
const vtkm::FloatDefault& Z,
const vtkm::FloatDefault& psi,
const XGCParams& Params)
{
constexpr vtkm::FloatDefault epsil_psi = 1e-5;
if ((psi <= (Params.eq_x_psi - epsil_psi) &&
-(R - Params.eq_x_r) * Params.eq_x_slope + (Z - Params.eq_x_z) > 0 &&
-(R - Params.eq_x2_r) * Params.eq_x2_slope + (Z - Params.eq_x2_z) < 0) ||
(psi > Params.eq_x_psi - epsil_psi && psi <= Params.eq_x2_psi - epsil_psi &&
-(R - Params.eq_x2_r) * Params.eq_x2_slope + (Z - Params.eq_x2_z) < 0) ||
psi > Params.eq_x2_psi - epsil_psi)
{
return true;
}
return false;
}
template <typename ScalarPortal>
VTKM_EXEC bool HighOrderB(const vtkm::Vec3f& pointRPZ,
ParticleMetaData& metadata,
const XGCParams& Params,
const ScalarPortal& Coeff_1D,
const ScalarPortal& Coeff_2D)
{
vtkm::FloatDefault R = pointRPZ[0], Z = pointRPZ[2];
vtkm::Id r_i = GetIndex(R, Params.nr, Params.eq_min_r, Params.dr_inv);
vtkm::Id z_i = GetIndex(Z, Params.nz, Params.zmin, Params.dz_inv);
vtkm::FloatDefault Rc = Params.eq_min_r + (vtkm::FloatDefault)(r_i)*Params.dr;
vtkm::FloatDefault Zc = Params.eq_min_z + (vtkm::FloatDefault)(z_i)*Params.dz;
auto Rc_1 = Rc + Params.dr;
auto Zc_1 = Zc + Params.dz;
Rc = (Rc + Rc_1) * 0.5;
Zc = (Zc + Zc_1) * 0.5;
//Get the coeffcients (z,r,4,4)
vtkm::Matrix<vtkm::FloatDefault, 4, 4> acoeff;
//offset = ri * nz + zi
//vtkm::Id offset = (r_i * ncoeff + z_i) * 16; //DRP
vtkm::Id offset = (z_i * Params.ncoeff_r + r_i) * 16;
vtkm::FloatDefault psi, dpsi_dr, dpsi_dz, d2psi_d2r, d2psi_drdz, d2psi_d2z;
EvalBicub2(
R, Z, Rc, Zc, offset, Coeff_2D, psi, dpsi_dr, dpsi_dz, d2psi_drdz, d2psi_d2r, d2psi_d2z);
if (psi < 0)
psi = 0;
metadata.Psi = psi;
//PSI = psi;
metadata.gradPsi_rzp[0] = dpsi_dr;
metadata.gradPsi_rzp[1] = dpsi_dz;
metadata.gradPsi_rzp[2] = 0;
metadata.dpsi_dr = dpsi_dr;
metadata.dpsi_dz = dpsi_dz;
metadata.d2psi_drdz = d2psi_drdz;
metadata.d2psi_d2r = d2psi_d2r;
metadata.d2psi_d2z = d2psi_d2z;
vtkm::FloatDefault fld_I, fld_dIdpsi;
bool isR12 = IsRegion12(R, Z, psi, Params);
if (!isR12)
{
fld_I = I_interpol(psi, 0, 3, Params, Coeff_1D);
fld_dIdpsi = I_interpol(psi, 1, 3, Params, Coeff_1D);
}
else
{
fld_I = I_interpol(psi, 0, 1, Params, Coeff_1D);
fld_dIdpsi = I_interpol(psi, 1, 1, Params, Coeff_1D);
}
vtkm::FloatDefault over_r = 1 / R;
vtkm::FloatDefault over_r2 = over_r * over_r;
vtkm::FloatDefault Br = -dpsi_dz * over_r;
vtkm::FloatDefault Bz = dpsi_dr * over_r;
vtkm::FloatDefault Bp = fld_I * over_r;
metadata.B0_rzp = vtkm::Vec3f(Br, Bz, Bp);
const vtkm::FloatDefault bp_sign = 1.0;
auto dBr_dr = (dpsi_dz * over_r2 - d2psi_drdz * over_r) * bp_sign;
auto dBr_dz = -d2psi_d2z * over_r * bp_sign;
auto dBr_dp = 0.0 * bp_sign;
auto dBz_dr = (-dpsi_dr * over_r2 + d2psi_d2r * over_r) * bp_sign;
auto dBz_dz = d2psi_drdz * over_r * bp_sign;
auto dBz_dp = 0.0 * bp_sign;
auto dBp_dr = dpsi_dr * fld_dIdpsi * over_r - fld_I * over_r2;
auto dBp_dz = fld_dIdpsi * dpsi_dz * over_r;
//auto dBp_dp = jacobian_rzp[2][2];
//calculate curl_B
//vtkm::Vec3f curlB_rzp;
metadata.curlB_rzp[0] = dBz_dp * over_r - dBp_dz;
metadata.curlB_rzp[1] = Bp * over_r + dBp_dr - dBr_dp * over_r;
metadata.curlB_rzp[2] = dBr_dz - dBz_dr;
//std::cout<<"curl_B_rzp= "<<curlB_rzp<<std::endl;
//calculate curl_nb
vtkm::FloatDefault Bmag = vtkm::Magnitude(metadata.B0_rzp);
vtkm::FloatDefault over_B = 1. / Bmag, over_B2 = over_B * over_B;
vtkm::FloatDefault dBdr, dBdz;
dBdr = (Br * dBr_dr + Bp * dBp_dr + Bz * dBz_dr) * over_B;
dBdz = (Br * dBr_dz + Bp * dBp_dz + Bz * dBz_dz) * over_B;
//vtkm::Vec3f curl_nb_rzp;
metadata.curl_nb_rzp[0] = metadata.curlB_rzp[0] * over_B + (Bp * dBdz) * over_B2;
metadata.curl_nb_rzp[1] = metadata.curlB_rzp[1] * over_B + (-Bp * dBdr) * over_B2;
metadata.curl_nb_rzp[2] = metadata.curlB_rzp[2] * over_B + (Bz * dBdr - Br * dBdz) * over_B2;
return true;
}
template <typename ScalarPortal>
VTKM_EXEC vtkm::Vec3f HighOrderBOnly(const vtkm::Vec3f& ptRPZ,
const XGCParams& Params,
const ScalarPortal& Coeff_1D,
const ScalarPortal& Coeff_2D)
{
vtkm::FloatDefault R = ptRPZ[0], Z = ptRPZ[2];
vtkm::Vec3f ptRZ(R, Z, 0);
int r_i = GetIndex(R, Params.nr, Params.rmin, Params.dr_inv);
int z_i = GetIndex(Z, Params.nz, Params.zmin, Params.dz_inv);
// rc(i), zc(j)
vtkm::FloatDefault Rc = Params.rmin + (vtkm::FloatDefault)(r_i)*Params.dr;
vtkm::FloatDefault Zc = Params.zmin + (vtkm::FloatDefault)(z_i)*Params.dz;
auto Rc_1 = Rc + Params.dr;
auto Zc_1 = Zc + Params.dz;
Rc = (Rc + Rc_1) * 0.5;
Zc = (Zc + Zc_1) * 0.5;
//Get the coeffcients (z,r,4,4)
vtkm::Matrix<vtkm::FloatDefault, 4, 4> acoeff;
//vtkm::Id offset = (r_i * ncoeff + z_i) * 16; //DRP
vtkm::Id offset = (z_i * Params.ncoeff_r + r_i) * 16;
/*
std::cout<<"InterpolatePsi: "<<vtkm::Vec2f(R,Z)<<std::endl;
std::cout<<" i/j= "<<r_i<<" "<<z_i<<std::endl;
std::cout<<" ncoeff= "<<ncoeff<<std::endl;
std::cout<<" offset= "<<offset<<std::endl;
std::cout<<" Rc/Zc= "<<Rc<<" "<<Zc<<std::endl;
*/
vtkm::FloatDefault psi, dpsi_dr, dpsi_dz, d2psi_d2r, d2psi_drdz, d2psi_d2z;
EvalBicub2(
R, Z, Rc, Zc, offset, Coeff_2D, psi, dpsi_dr, dpsi_dz, d2psi_drdz, d2psi_d2r, d2psi_d2z);
if (psi < 0)
psi = 0;
bool isR12 = IsRegion12(R, Z, psi, Params);
vtkm::FloatDefault fld_I;
if (!isR12)
fld_I = I_interpol(psi, 0, 3, Params, Coeff_1D);
else
fld_I = I_interpol(psi, 0, 0, Params, Coeff_1D);
vtkm::FloatDefault over_r = 1 / R;
vtkm::FloatDefault Br = -dpsi_dz * over_r;
vtkm::FloatDefault Bz = dpsi_dr * over_r;
vtkm::FloatDefault Bp = fld_I * over_r;
return vtkm::Vec3f(Br, Bz, Bp);
}
template <typename ScalarPortal>
VTKM_EXEC vtkm::Vec3f GetB(vtkm::Vec3f& pt_rpz,
const XGCParams& Params,
const ScalarPortal& Coeff_1D,
const ScalarPortal& Coeff_2D)
{
auto B0_rzp = HighOrderBOnly(pt_rpz, Params, Coeff_1D, Coeff_2D);
vtkm::Vec3f B0_rpz(B0_rzp[0], B0_rzp[2], B0_rzp[1]);
// This gives is the time derivative: Br = dR/dt, Bz= dZ/dt, B_phi/R = dphi/dt
// We need with respect to phi:
// dR/dphi = dR/dt / (dphi/dt) = Br / B_phi * R
// same for z;
B0_rpz[0] /= B0_rzp[2];
B0_rpz[2] /= B0_rzp[2];
B0_rpz[0] *= pt_rpz[0];
B0_rpz[2] *= pt_rpz[0];
return B0_rpz;
}
template <typename ScalarPortal>
VTKM_EXEC bool CalcFieldFollowingPt(const vtkm::Vec3f& pt_rpz,
const vtkm::Vec3f& B0_rpz,
const vtkm::FloatDefault& Phi0,
const vtkm::FloatDefault& Phi1,
const XGCParams& Params,
const ScalarPortal& coeff_1D,
const ScalarPortal& coeff_2D,
vtkm::Vec3f& x_ff_rpz)
{
using Ray3f = vtkm::Ray<vtkm::FloatDefault, 3, true>;
vtkm::FloatDefault R = pt_rpz[0];
vtkm::FloatDefault Phi = pt_rpz[1];
vtkm::FloatDefault Z = pt_rpz[2];
vtkm::FloatDefault PhiMid = Phi0 + (Phi1 - Phi0) / 2.0;
vtkm::Plane<> midPlane({ 0, PhiMid, 0 }, { 0, 1, 0 });
Ray3f ray_rpz({ R, Phi, Z }, B0_rpz);
//Get point on mid plane. Use the R,Z for this point for triangle finds.
vtkm::FloatDefault RP_T;
vtkm::Vec3f ptOnMidPlane_rpz;
bool b;
midPlane.Intersect(ray_rpz, RP_T, ptOnMidPlane_rpz, b);
//Now, do it using RK4 and two steps.
vtkm::FloatDefault h = (PhiMid - Phi) / 2.0;
//h = -h;
vtkm::FloatDefault h_2 = h / 2.0;
//k1 = F(p)
//k2 = F(p+hk1/2)
//k3 = F(p+hk2/2)
//k4 = F(p+hk3)
//Yn+1 = Yn + 1/6 h (k1+2k2+2k3+k4)
vtkm::Vec3f p0 = { R, Phi, Z }; //pt_rpz;
vtkm::Vec3f tmp, k1, k2, k3, k4;
//std::cout<<" p0 = "<<p0<<std::endl;
for (int i = 0; i < 2; i++)
{
k1 = GetB(p0, Params, coeff_1D, coeff_2D);
tmp = p0 + k1 * h_2;
k2 = GetB(tmp, Params, coeff_1D, coeff_2D);
tmp = p0 + k2 * h_2;
k3 = GetB(tmp, Params, coeff_1D, coeff_2D);
tmp = p0 + k3 * h;
k4 = GetB(tmp, Params, coeff_1D, coeff_2D);
vtkm::Vec3f vec = (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;
p0 = p0 + h * vec;
}
x_ff_rpz = p0;
return true;
}
template <typename PortalType>
VTKM_EXEC typename PortalType::ValueType Eval(const PortalType& portal,
const vtkm::Id& offset,
const vtkm::Vec3f& param,
const vtkm::Vec<vtkm::Id, 3>& vId)
{
#if 1
//Hard code...
const auto& v0 = portal.Get(vId[0] + offset);
const auto& v1 = portal.Get(vId[1] + offset);
const auto& v2 = portal.Get(vId[2] + offset);
const auto& w1 = param[0];
const auto& w2 = param[1];
const auto w0 = 1 - (w1 + w2);
auto ret = v0 * w0 + v1 * w1 + v2 * w2;
#else
T ret;
vtkm::VecVariable<T, 3> vals;
vals.Append(portal.Get(vId[0] + offset));
vals.Append(portal.Get(vId[1] + offset));
vals.Append(portal.Get(vId[2] + offset));
vtkm::exec::CellInterpolate(vals, param, vtkm::CellShapeTagTriangle(), ret);
#endif
return ret;
}
VTKM_EXEC inline void GetPlaneIdx(const vtkm::FloatDefault& phi,
const vtkm::FloatDefault& Period,
const XGCParams& Params,
vtkm::FloatDefault& phiN,
vtkm::Id& plane0,
vtkm::Id& plane1,
vtkm::FloatDefault& phi0,
vtkm::FloatDefault& phi1,
vtkm::Id& numRevs,
vtkm::FloatDefault& T)
{
numRevs = vtkm::Floor(vtkm::Abs(phi / Period));
phiN = phi;
if (phi < 0)
{
phiN += ((1 + numRevs) * Period);
}
else if (phi > Period)
{
phiN -= (numRevs * Period);
}
plane0 = vtkm::Floor(phiN / Params.dPhi);
if (plane0 == Params.NumPlanes)
plane0 = 0;
plane1 = plane0 + 1;
phi0 = static_cast<vtkm::FloatDefault>(plane0) * Params.dPhi;
phi1 = static_cast<vtkm::FloatDefault>(plane1) * Params.dPhi;
if (plane1 == Params.NumPlanes)
plane1 = 0;
T = (phiN - phi0) / (phi1 - phi0);
}
}
}
} //vtkm::worklet::flow
#endif //vtkm_filter_flow_worklet_XGCHelper_h

281
vtkm/filter/flow/worklet/XGCPoincare.h
View File

@ -1,281 +0,0 @@
//=============================================================================
//
// 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 vtkm_worklet_particleadvection_xgcpoincare_h
#define vtkm_worklet_particleadvection_xgcpoincare_h
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/filter/flow/worklet/XGCHelper.h>
namespace vtkm
{
namespace worklet
{
namespace flow
{
template <typename ParticleType>
class XGCPoincareExec
{
public:
VTKM_EXEC_CONT
XGCPoincareExec()
: NumParticles(0)
, MaxSteps(0)
, MaxPunctures(0)
, PunctureCounts()
, Validity()
, PunctureIDs()
, OutputR()
, OutputZ()
, OutputTheta()
, OutputPsi()
{
}
VTKM_CONT
XGCPoincareExec(vtkm::Id numParticles,
vtkm::Id maxSteps,
vtkm::Id maxPunctures,
const XGCParams& params,
vtkm::FloatDefault period,
const vtkm::cont::ArrayHandle<vtkm::FloatDefault>& coeff_1D,
const vtkm::cont::ArrayHandle<vtkm::FloatDefault>& coeff_2D,
const vtkm::cont::ArrayHandle<vtkm::Id>& punctureCounts,
const vtkm::cont::ArrayHandle<vtkm::Id>& validity,
const vtkm::cont::ArrayHandle<vtkm::Id>& punctureIDs,
const vtkm::cont::ArrayHandle<vtkm::FloatDefault>& outputR,
const vtkm::cont::ArrayHandle<vtkm::FloatDefault>& outputZ,
const vtkm::cont::ArrayHandle<vtkm::FloatDefault>& outputTheta,
const vtkm::cont::ArrayHandle<vtkm::FloatDefault>& outputPsi,
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token)
: NumParticles(numParticles)
, MaxSteps(maxSteps)
, MaxPunctures(maxPunctures)
, Params(params)
, Period(period)
{
Coeff_1D = coeff_1D.PrepareForInput(device, token);
Coeff_2D = coeff_2D.PrepareForInput(device, token);
PunctureCounts = punctureCounts.PrepareForInPlace(device, token);
Validity = validity.PrepareForInPlace(device, token);
PunctureIDs = punctureIDs.PrepareForInPlace(device, token);
OutputR = outputR.PrepareForOutput(this->NumParticles * this->MaxPunctures, device, token);
OutputZ = outputZ.PrepareForOutput(this->NumParticles * this->MaxPunctures, device, token);
OutputTheta =
outputTheta.PrepareForOutput(this->NumParticles * this->MaxPunctures, device, token);
OutputPsi = outputPsi.PrepareForOutput(this->NumParticles * this->MaxPunctures, device, token);
}
//template <typename ParticleType>
VTKM_EXEC void Analyze(const vtkm::Id index,
const ParticleType& oldParticle,
ParticleType& newParticle) const
{
vtkm::Id numRevs0 = vtkm::Floor(vtkm::Abs(oldParticle.Pos[1] / this->Period));
vtkm::Id numRevs1 = vtkm::Floor(vtkm::Abs(newParticle.Pos[1] / this->Period));
//if (this->SaveTraces)
// Intersections.Set(idx*this->MaxIter + particle.NumSteps, particle.Pos);
//std::cout<<std::setprecision(12)<<" Step: "<<particle.Pos<<std::endl;
if (numRevs1 > numRevs0)
{
auto R = newParticle.Pos[0];
auto Z = newParticle.Pos[2];
auto theta = vtkm::ATan2(Z - this->Params.eq_axis_z, R - this->Params.eq_axis_r);
if (theta < 0)
theta += vtkm::TwoPi();
//calcualte psi. need to stash psi on the particle somehow....
vtkm::Vec3f ptRPZ = newParticle.Pos;
ParticleMetaData metadata;
HighOrderB(ptRPZ, metadata, this->Params, this->Coeff_1D, this->Coeff_2D);
vtkm::Id loc = (index * this->MaxPunctures) + newParticle.NumPunctures;
this->Validity.Set(loc, 1);
this->OutputR.Set(loc, R);
this->OutputZ.Set(loc, Z);
this->OutputTheta.Set(loc, theta);
this->OutputPsi.Set(loc, metadata.Psi / this->Params.eq_x_psi);
this->PunctureIDs.Set(loc, index);
newParticle.NumPunctures++;
}
}
private:
using ScalarPortalR = typename vtkm::cont::ArrayHandle<vtkm::FloatDefault>::ReadPortalType;
using ScalarPortalW = typename vtkm::cont::ArrayHandle<vtkm::FloatDefault>::WritePortalType;
using IdPortal = typename vtkm::cont::ArrayHandle<vtkm::Id>::WritePortalType;
using VecPortal = typename vtkm::cont::ArrayHandle<vtkm::Vec3f>::WritePortalType;
vtkm::Id NumParticles;
vtkm::Id MaxSteps;
vtkm::Id MaxPunctures;
vtkm::FloatDefault Period;
XGCParams Params;
ScalarPortalR Coeff_1D;
ScalarPortalR Coeff_2D;
ScalarPortalW OutputR;
ScalarPortalW OutputZ;
ScalarPortalW OutputTheta;
ScalarPortalW OutputPsi;
IdPortal PunctureIDs;
//VecPortal Intersections;
IdPortal PunctureCounts;
IdPortal Validity;
};
template <typename ParticleType>
class XGCPoincare : public vtkm::cont::ExecutionObjectBase
{
public:
// Intended to store advected particles after Finalize
vtkm::cont::ArrayHandle<ParticleType> Particles;
vtkm::cont::ArrayHandle<vtkm::Id> PunctureIDs;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> OutputR;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> OutputZ;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> OutputTheta;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> OutputPsi;
//Helper functor for compacting history
struct IsOne
{
template <typename T>
VTKM_EXEC_CONT bool operator()(const T& x) const
{
return x == T(1);
}
};
VTKM_CONT
XGCPoincare(vtkm::Id maxSteps,
vtkm::Id maxPunctures,
const vtkm::cont::ArrayHandle<vtkm::FloatDefault>& coeff_1D,
const vtkm::cont::ArrayHandle<vtkm::FloatDefault>& coeff_2D,
const XGCParams& params,
vtkm::FloatDefault period)
: Particles()
, MaxSteps(maxSteps)
, MaxPunctures(maxPunctures)
, Coeff_1D(coeff_1D)
, Coeff_2D(coeff_2D)
, Params(params)
, Period(period)
{
}
VTKM_CONT
void UseAsTemplate(const XGCPoincare& other)
{
this->MaxSteps = other.MaxSteps;
this->MaxPunctures = other.MaxPunctures;
this->Period = other.Period;
this->Params = other.Params;
this->Coeff_1D = other.Coeff_1D;
this->Coeff_2D = other.Coeff_2D;
}
VTKM_CONT XGCPoincareExec<ParticleType> PrepareForExecution(vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
{
return XGCPoincareExec<ParticleType>(this->NumParticles,
this->MaxSteps,
this->MaxPunctures,
this->Params,
this->Period,
this->Coeff_1D,
this->Coeff_2D,
this->PunctureCounts,
this->Validity,
this->PunctureIDs,
this->OutputR,
this->OutputZ,
this->OutputTheta,
this->OutputPsi,
device,
token);
}
VTKM_CONT bool SupportPushOutOfBounds() const { return false; }
VTKM_CONT
void InitializeAnalysis(const vtkm::cont::ArrayHandle<ParticleType>& particles)
{
this->NumParticles = particles.GetNumberOfValues();
//Create Validity initialized to zero.
vtkm::cont::ArrayHandleConstant<vtkm::Id> validity(0, this->NumParticles * this->MaxPunctures);
vtkm::cont::ArrayCopy(validity, this->Validity);
//Create PunctureCount initialized to zero.
vtkm::cont::ArrayHandleConstant<vtkm::Id> punctureCounts(0, this->NumParticles);
vtkm::cont::ArrayCopy(punctureCounts, this->PunctureCounts);
vtkm::cont::ArrayHandleConstant<vtkm::Id> initIds(-1, this->NumParticles * this->MaxPunctures);
vtkm::cont::ArrayCopy(initIds, this->PunctureIDs);
}
VTKM_CONT
void FinalizeAnalysis(vtkm::cont::ArrayHandle<ParticleType>& particles)
{
(void)particles;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> outR, outZ, outTheta, outPsi;
vtkm::cont::ArrayHandle<vtkm::Id> outIds;
vtkm::cont::Algorithm::CopyIf(this->OutputR, this->Validity, outR, IsOne());
vtkm::cont::Algorithm::CopyIf(this->OutputZ, this->Validity, outZ, IsOne());
vtkm::cont::Algorithm::CopyIf(this->OutputTheta, this->Validity, outTheta, IsOne());
vtkm::cont::Algorithm::CopyIf(this->OutputPsi, this->Validity, outPsi, IsOne());
vtkm::cont::Algorithm::CopyIf(this->PunctureIDs, this->Validity, outIds, IsOne());
vtkm::cont::ArrayCopyShallowIfPossible(outR, this->OutputR);
vtkm::cont::ArrayCopyShallowIfPossible(outZ, this->OutputZ);
vtkm::cont::ArrayCopyShallowIfPossible(outTheta, this->OutputTheta);
vtkm::cont::ArrayCopyShallowIfPossible(outPsi, this->OutputPsi);
vtkm::cont::ArrayCopyShallowIfPossible(outIds, this->PunctureIDs);
}
VTKM_CONT static bool MakeDataSet(vtkm::cont::DataSet& dataset,
const std::vector<XGCPoincare>& results)
{
(void)dataset;
(void)results;
return true;
}
private:
vtkm::Id NumParticles;
vtkm::Id MaxSteps;
vtkm::Id MaxPunctures;
// For Higher order interpolation
// on RZ uniform grid, not the triangle grid.
vtkm::cont::ArrayHandle<vtkm::FloatDefault> Coeff_1D;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> Coeff_2D;
XGCParams Params;
vtkm::FloatDefault Period;
// Output Arrays
vtkm::cont::ArrayHandle<vtkm::Id> PunctureCounts;
//vtkm::cont::ArrayHandle<vtkm::Id> InitialPunctureCounts;
vtkm::cont::ArrayHandle<vtkm::Id> Validity;
};
} // namespace particleadvection
} // namespace worklet
} // namespace vtkm
#endif //vtkm_worklet_particleadvection_xgcpoincare_h