1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
|
/* Copyright 2016-2020 Maxence Thevenet, Yinjian Zhao
*
* This file is part of WarpX.
*
* License: BSD-3-Clause-LBNL
*/
#ifndef FILTERFUNCTORS_H
#define FILTERFUNCTORS_H
#include "Particles/Pusher/GetAndSetPosition.H"
#include "Particles/WarpXParticleContainer.H"
#include "Utils/WarpXConst.H"
#include <AMReX_Gpu.H>
#include <AMReX_Parser.H>
#include <AMReX_Random.H>
using SuperParticleType = typename WarpXParticleContainer::SuperParticleType;
/**
* \brief Used to keep track of what inputs units a filter function should expect.
* "WarpX units" means the momentum is "gamma*v" (aka proper velocity)
* "SI" means the momentum is mass*gamma*v.
*/
enum struct InputUnits{WarpX, SI};
/**
* \brief Functor that returns 0 or 1 depending on a random draw per particle
*/
struct RandomFilter
{
/** constructor
* \param a_is_active whether the test is active
* \param a_fraction fraction of particles to select
*/
RandomFilter(bool a_is_active, amrex::Real a_fraction)
: m_is_active(a_is_active), m_fraction(a_fraction) {}
/**
* \brief draw random number, return 1 if number < m_fraction, 1 otherwise
* \param p one particle
* \param engine the random number state and factory
* \return whether or not the particle is selected
*/
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
bool operator () (const SuperParticleType& p, const amrex::RandomEngine& engine) const noexcept
{
amrex::ignore_unused(p);
return ( (!m_is_active) || (amrex::Random(engine) < m_fraction) );
}
private:
const bool m_is_active; //! select all particles if false
const amrex::Real m_fraction = 1.0; //! range: [0.0:1.0] where 0 is no & 1 is all particles
};
/**
* \brief Functor that returns 1 if stride divide particle_id, 0 otherwise
*/
struct UniformFilter
{
/** constructor
* \param a_is_active whether the test is active
* \param a_stride one particle every a_stride is written to file
*/
UniformFilter(bool a_is_active, int a_stride)
: m_is_active(a_is_active), m_stride(a_stride) {}
/**
* \brief return 1 if stride divide particle_id, 0 otherwise
* \param p one particle
* \return whether or not the particle is selected
*/
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
bool operator () (const SuperParticleType& p, const amrex::RandomEngine&) const noexcept
{
return ( (!m_is_active) || ( p.id()%m_stride == 0 ) );
}
private:
const bool m_is_active; //! select all particles if false
const int m_stride = 0; //! selection of every n-th particle
};
/**
* \brief Functor that returns 0 or 1 depending on a parser selection
*/
struct ParserFilter
{
/** constructor
* \param a_is_active whether the test is active
* \param a_filter_parser parser taking t, x, y, z, ux, uy, and uz, and returning a boolean for selected particle
* \param a_mass mass of the particle species
* \param time simulation time on the coarsest level
*/
ParserFilter(bool a_is_active, amrex::ParserExecutor<7> const& a_filter_parser,
const amrex::ParticleReal a_mass, const amrex::Real time):
m_is_active{a_is_active},
m_function_partparser{a_filter_parser},
m_mass{a_mass},
m_t{time},
m_units{InputUnits::WarpX}{}
/**
* \brief return 1 if the particle is selected by the parser
* \param p one particle
* \return whether or not the particle is selected
*/
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
bool operator () (const SuperParticleType& p, const amrex::RandomEngine&) const noexcept
{
using namespace amrex::literals;
if (!m_is_active){
return true;
}
else{
amrex::ParticleReal x, y, z;
get_particle_position(p, x, y, z);
amrex::Real ux = p.rdata(PIdx::ux)/PhysConst::c;
amrex::Real uy = p.rdata(PIdx::uy)/PhysConst::c;
amrex::Real uz = p.rdata(PIdx::uz)/PhysConst::c;
if (m_units == InputUnits::SI)
{
ux /= m_mass;
uy /= m_mass;
uz /= m_mass;
}
// ux, uy, uz are now in beta*gamma
// This is actually a binary true/false (1/0) check,
// but the parser returns floating point types
return (m_function_partparser(m_t,x,y,z,ux,uy,uz) != 0.0_rt);
}
}
private:
/** Whether this diagnostics is activated. Select all particles if false */
const bool m_is_active;
public:
/** Parser function with 7 input variables, t,x,y,z,ux,uy,uz */
amrex::ParserExecutor<7> const m_function_partparser;
/** Mass of particle species */
amrex::ParticleReal m_mass;
/** Store physical time on the coarsest level. */
amrex::Real m_t;
/** keep track of momentum units particles will come in with **/
InputUnits m_units;
};
/**
* \brief Functor that returns 1 if the particle is inside a given axis-aligned region
* defined by amrex::RealBox, 0 otherwise.
*/
struct GeometryFilter
{
GeometryFilter(bool a_is_active, amrex::RealBox a_domain)
: m_is_active(a_is_active), m_domain(a_domain) {}
/**
* \brief return 1 if the partcile is within the region described by the RealBox
* \param p one particle
* \return whether or not the particle is inside the region defined by m_domain
*/
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
bool operator () (const SuperParticleType& p, const amrex::RandomEngine&) const noexcept
{
if ( !m_is_active ) return 1;
return ! (AMREX_D_TERM( (p.pos(0) < m_domain.lo(0)) || (p.pos(0) > m_domain.hi(0) ),
|| (p.pos(1) < m_domain.lo(1)) || (p.pos(1) > m_domain.hi(1) ),
|| (p.pos(2) < m_domain.lo(2)) || (p.pos(2) > m_domain.hi(2) )));
}
private:
/** Whether this diagnostics is activated. Select all particles if false */
const bool m_is_active;
/** Physical extent of the axis-aligned region used for particle check */
const amrex::RealBox m_domain;
};
#endif // FILTERFUNCTORS_H
|
