Particle Filters

Section author: Axel Huebl, Sergei Bastrakov

A common task in both modeling, initializing and in situ processing (output) is the selection of particles of a particle species by attributes. PIConGPU implements such selections as particle filters.

Particle filters are simple mappings assigning each particle of a species either true or false (ignore / filter out). These filters can be defined in particleFilters.param.

Example

Let us select particles with momentum vector within a cone with an opening angle of five degrees (pinhole):

namespace picongpu
{
namespace particles
{
namespace filter
{
    struct FunctorParticlesForwardPinhole
    {
        static constexpr char const * name = "forwardPinhole";

        template< typename T_Particle >
        HDINLINE bool operator()(
            T_Particle const & particle
        )
        {
            bool result = false;
            float3_X const mom = particle[ momentum_ ];
            float_X const normMom = pmacc::math::l2norm( mom );

            if( normMom > float_X( 0. ) )
            {
                /* place detector in y direction, "infinite distance" to target,
                 * and five degree opening angle
                 */
                constexpr float_X openingAngle = 5.0 * PI / 180.;
                float_X const dotP = mom.y() / normMom;
                float_X const degForw = math::acos( dotP );

                if( math::abs( degForw ) <= openingAngle * float_X( 0.5 ) )
                    result = true;
            }
            return result;
        }
    };
    using ParticlesForwardPinhole = generic::Free<
       FunctorParticlesForwardPinhole
    >;

and add ParticlesForwardPinhole to the AllParticleFilters list:

    using AllParticleFilters = MakeSeq_t<
        All,
        ParticlesForwardPinhole
    >;

} // namespace filter
} // namespace particles
} // namespace picongpu

Filtering Based on Global Position

A particle only stores its relative position inside the cell. Thus, with a simple filter like one shown above, it is not possible to get global position of a particle. However, there are helper wrapper filters that provide such information in addition to the particle data.

For a special case of slicing along one axis this is a simple existing filter that only needs to be parametrized:

namespace picongpu
{
namespace particles
{
namespace filter
{
    namespace detail
    {
        //! Parameters to be used with RelativeGlobalDomainPosition, change the values inside
        struct SliceParam
        {
            // Lower bound in relative coordinates: global domain is [0.0, 1.0]
            static constexpr float_X lowerBound = 0.55_X;

            // Upper bound in relative coordinates
            static constexpr float_X upperBound = 0.6_X;

            // Axis: x = 0; y= 1; z = 2
            static constexpr uint32_t dimension = 0;

            // Text name of the filter, will be used in .cfg file
            static constexpr char const* name = "slice";
        };

        //! Use the existing RelativeGlobalDomainPosition filter with our parameters
        using Slice = RelativeGlobalDomainPosition<SliceParam>;
    }

and add detail::Slice to the AllParticleFilters list:

    using AllParticleFilters = MakeSeq_t<
        All,
        detail::Slice
    >;

} // namespace filter
} // namespace particles
} // namespace picongpu

For a more general case of filtering based on cell index (possibly combined with other particle properties) use the following pattern:

namespace picongpu
{
namespace particles
{
namespace filter
{
    namespace detail
    {
        struct AreaFilter
        {
            static constexpr char const* name = "areaFilter";

            template<typename T_Particle>
            HDINLINE bool operator()(
                DataSpace<simDim> const totalCellOffset,
                T_Particle const & particle
            )
            {
                /* Here totalCellOffset is the cell index of the particle in the total coordinate system.
                 * So we can define conditions based on both cell index and other particle data.
                 */
                return (totalCellOffset.x() >= 10) && (particle[momentum_].x() < 0.0_X);
             }
         };

         //! Wrap AreaFilter so that it fits the general filter interface
         using Area = generic::FreeTotalCellOffset<AreaFilter>;
    }

and add detail::Area to the AllParticleFilters list:

    using AllParticleFilters = MakeSeq_t<
        All,
        detail::Area
    >;

} // namespace filter
} // namespace particles
} // namespace picongpu

Limiting Filters to Eligible Species

Besides the list of pre-defined filters with parametrization, users can also define generic, “free” implementations as shown above. All filters are added to AllParticleFilters and then combined with all available species from VectorAllSpecies (see speciesDefinition.param).

In the case of user-defined free filters we can now check if each species in VectorAllSpecies fulfills the requirements of the filter. That means: if one accesses specific attributes or flags of a species in a filter, they must exist or will lead to a compile error.

As an example, probe particles usually do not need a momentum attribute which would be used for an energy filter. So they should be ignored from compilation when combining filters with particle species.

In order to exclude all species that have no momentum attribute from the ParticlesForwardPinhole filter, specialize the C++ trait SpeciesEligibleForSolver. This trait is implemented to be checked during compile time when combining filters with species:

// ...

} // namespace filter

namespace traits
{
    template<
        typename T_Species
    >
    struct SpeciesEligibleForSolver<
        T_Species,
        filter::ParticlesForwardPinhole
    >
    {
        using type = typename pmacc::traits::HasIdentifiers<
            typename T_Species::FrameType,
            MakeSeq_t< momentum >
        >::type;
    };
} // namespace traits
} // namespace particles
} // namespace picongpu