XGC1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
update_particle_flags.hpp
Go to the documentation of this file.
1 #ifndef UPDATE_PARTICLE_FLAGS_HPP
2 #define UPDATE_PARTICLE_FLAGS_HPP
3 
4 #include "particle_stream.hpp"
5 #include "magnetic_field.hpp"
6 #include "push_controls.hpp"
7 
8 #ifdef ESC_PTL
9 // In the following subroutine, following flags need to be set.
10 // - is_escaped
11 // - is_divertor --> after sheath
12 // - is_outboard
13 // - is_to_write : when is_escpaed or is_divertor set
14 // - was_inside : set after it is used to check is_escaped
15 // - do not consider the situation that is_escpaed and is_divertor happen together.
16 template<class Device>
17 KOKKOS_INLINE_FUNCTION void update_particle_flag1(const PushControls &push_controls, const MagneticField<Device> &magnetic_field, SimdParticles &part_tmp, Simd<double> &psi)
18 {
20  Simd<bool> bt0, bt1, bt2, bt3; // temp. boolean
21 
22  for (int i_simd = 0; i_simd < SIMD_SIZE; i_simd++){
23  // zero out above 25 bits. (8 bits for flags and 17 bits for esc_step)
24  part_tmp.flag[i_simd] = part_tmp.flag[i_simd] & ps_flag_step_mask;
25 
26  bt0[i_simd] = flag_check(part_tmp.flag[i_simd], ps_is_to_write); // it needs to be written (determined previously), do not change any flag
27  //bt1[i_simd] = psi[i_simd] < XPOINT_PSI && x.z[i_simd] > XPOINT_Z ; // new position is inside
28  bt1[i_simd] = (magnetic_field.is_in_region_1(part_tmp.ph.r[i_simd],part_tmp.ph.z[i_simd],psi[i_simd])); // new position is in region 1 ( = inside)
29  bt2[i_simd] = (!bt1[i_simd]) && flag_check(part_tmp.flag[i_simd], ps_was_inside); // is_outside && was_inside --> escaped
30  bt3[i_simd] = part_tmp.ph.r[i_simd] > magnetic_field.equil.axis_r ; // is_outboard
31 
32  // set is_escpaed, is_outboard, is_to_write (only when bt0 false)
33  itmp[i_simd] = part_tmp.flag[i_simd];
34  flag_set(itmp[i_simd], bt2[i_simd], ps_is_escaped); // set is_escpaed as bt2
35  flag_set(itmp[i_simd], bt2[i_simd], ps_is_to_write); // set is_to_write (always when is_escaped true.)
36  flag_set(itmp[i_simd], bt3[i_simd], ps_is_outboard); // set is_outboard as bt3
37  flag_set(itmp[i_simd], bt1[i_simd], ps_was_inside); // set was_insde to future poistion
38 
39  // Particles contributing to diffusion coefficient analysis are marked with ps_is_tracked.
40  // Only particles in the scrape-off layer contribute to this calculation.
41  // First, particles crossing the separatrix into the scrape-off layer need to have the ps_was_traced flag set
42  // Condition: ps_is_escaped == true (set above)
43  flag_set(itmp[i_simd], bt2[i_simd] || flag_check(itmp[i_simd],ps_is_tracked), ps_is_tracked);
44 
45  }
46 
47  //Set escaped time step -- separate loop due to if stetement. -- SIMD might not work.
48  for (int i_simd = 0; i_simd < SIMD_SIZE; i_simd++){
49  if(bt2[i_simd]) {
50  //escaped flag set true
51  long long int l_step = push_controls.gstep; // gstep in long long int data type.
52  itmp[i_simd] = itmp[i_simd] & ps_flag_mask; // zero out above 8 bits. Clearing old value.
53  itmp[i_simd] += (l_step << ps_flag_bits); // makes 8 zeros and adding to itmp.
54  // if l_step is 131072 (2^17) or above, overflow happens. Currently XGC output supports up to 99,999 time step in output name.
55  }
56  }
57 
58  // update particle flag only when is_to_write == false
59  for (int i_simd = 0; i_simd < SIMD_SIZE; i_simd++){
60  part_tmp.flag[i_simd] = (bt0[i_simd]) ? (part_tmp.flag[i_simd]) : (itmp[i_simd]);
61  }
62 
63  return;
64 }
65 
66 // For Gordon Bell 2024, the condition for tracked particles is changed.
67 // When a particle is inside of a flux shell (finite depth like 0.8 < psi < 0.81) at a certain time step,
68 // is_to_write and is_tracked will be true. Otherwise, it will be false.
69 // This update will not be called for a several time steps.
70 // It only needs is_to_write
71 template<class Device>
72 KOKKOS_INLINE_FUNCTION void update_particle_flag_diff(const PushControls &push_controls, const MagneticField<Device> &magnetic_field, SimdParticles &part_tmp, Simd<double> &psi, int idx, int i_cycle)
73 {
75  Simd<bool> bt0; // temp. boolean
76 
77  if(push_controls.ipc==2 && push_controls.gstep%push_controls.particle_stream.reset_period==0){
78  for (int i_simd = 0; i_simd < SIMD_SIZE; i_simd++){
79  // zero out above 25 bits. (8 bits for flags and 17 bits for esc_step)
80  part_tmp.flag[i_simd] = part_tmp.flag[i_simd] & ps_flag_step_mask;
81 
82 
83  bt0[i_simd] = (push_controls.particle_stream.pin < psi[i_simd]) && ( psi[i_simd]<push_controls.particle_stream.pout);
84 
85  // set is_to_write
86  itmp[i_simd] = part_tmp.flag[i_simd];
87  flag_set(itmp[i_simd], bt0[i_simd], ps_is_to_write); // set is_to_write
88  flag_set(itmp[i_simd], bt0[i_simd], ps_is_tracked); // set is_to_write
89 
90  // update particle flag
91  part_tmp.flag[i_simd] = itmp[i_simd];
92 
93  }
94  }
95 
96  return;
97 }
98 
99 KOKKOS_INLINE_FUNCTION void update_particle_flag2(Simd<bool> &not_in_triangle, Simd<double> &dw, SimdParticles &part_tmp)
100 {
101  Simd<long long int> itmp;
102  Simd<bool> bt0; // temp. boolean
103 
104  // update divertor flag (+ is_to_write) when the particle hit the divertor
105  // only send the particle with ipc==2 loop.
106  // can this be better to move after sheath_calculaton??
107  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
108 
109  bt0[i_simd] = flag_check(part_tmp.flag[i_simd], ps_is_to_write); // it needs to be written (determined previously), do not change any flag
110 
111  //beware that the following is called only when at_least_one is true
112  bool dw_nonzero = not_in_triangle[i_simd] & (dw[i_simd]!=0.0);
113 
114  itmp[i_simd]=part_tmp.flag[i_simd];
115  flag_set(itmp[i_simd], dw_nonzero, ps_is_divertor); // set is_divertor as dw_nonzero
116  flag_set(itmp[i_simd], dw_nonzero, ps_is_to_write); // set is_to_write the same as is_to_write
117 
118  part_tmp.flag[i_simd] = (bt0[i_simd] || !dw_nonzero) ? (part_tmp.flag[i_simd]) : (itmp[i_simd]); // update it only is_to_write == false & dw_nonzero == true
119  }
120 
121  // Pack dw value to flag
122  //beware that the following is called only when at_least_one is true
123  // flag has zero bits above flag_mask (3 bytes + 1bit) at this stage.
124  // dw can use 38 bits. dw * 2^nb1 will be saved as integer value.
125  // Let dw be between min_dw and max_dw
126  // This should not exceed 2^(63-nb1-nb2)
127  constexpr double min_dw=-256.;
128  constexpr double max_dw= 256.;
129 
130  // following factor1 and factor2 should be consistent unpacking fortran routine in
131  constexpr int nb1=ps_dw_factor_bits; // nb1 is mutiplication factor to make floating point to integer
132  constexpr int nb2=ps_flag_step_bits; // nb2 is mutiplication factor to clear last 8 + 17 bits.
133  constexpr int factor1 = 1u << nb1; // 2^28 , not more than 31
134  constexpr int factor2= 1u << nb2; // 2^25
135 
136  // loop over all particles to utilize SIMD
137  // dw should be zero for the particles that did not hit the divertor
138  // zero is added at the end.
139  // even though it is not zero, adding dw should be harmless.
140  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
141  dw[i_simd] = min(max_dw, max(min_dw, dw[i_simd]));
142  itmp[i_simd] = dw[i_simd] * static_cast<double>(factor1); // itmp is the integer to place at first 38 bits.
143  itmp[i_simd] = itmp[i_simd] * factor2; // multiplication is used to conserve sign, instead of '<< nb2'
144 
145  // 2022-03-03: Removed, RK
146  // 2022-03-21: Revived.
147  part_tmp.flag[i_simd]+=itmp[i_simd]; // add dw and flag. itmp should have zeros in the last 8 bit.
148  }
149  return;
150 }
151 #endif
152 #endif
constexpr int ps_flag_bits
Definition: particle_stream.hpp:49
KOKKOS_INLINE_FUNCTION bool is_in_region_1(double r, double z, double psi) const
Definition: magnetic_field.tpp:73
constexpr long long int ps_is_escaped
Definition: particle_stream.hpp:39
constexpr int ps_dw_factor_bits
Definition: particle_stream.hpp:52
constexpr int ps_flag_step_bits
Definition: particle_stream.hpp:51
KOKKOS_INLINE_FUNCTION void flag_set(long long int &flags, const bool tf, const long long int flag_loc)
Definition: particle_stream.tpp:9
constexpr long long int ps_is_to_write
Definition: particle_stream.hpp:36
Definition: push_controls.hpp:9
Definition: magnetic_field.hpp:13
Equilibrium equil
The object containing information about the magnetic equilibrium.
Definition: magnetic_field.hpp:45
constexpr long long int ps_flag_mask
Definition: particle_stream.hpp:55
KOKKOS_INLINE_FUNCTION bool flag_check(const long long int flags, const long long int flag_loc)
Definition: particle_stream.tpp:15
constexpr long long int ps_is_outboard
Definition: particle_stream.hpp:42
Simd< double > r
Definition: particles.hpp:18
constexpr long long int ps_is_divertor
Definition: particle_stream.hpp:41
constexpr long long int ps_flag_step_mask
Definition: particle_stream.hpp:54
SimdPhase ph
Definition: particles.hpp:59
Definition: particles.hpp:58
int gstep
Definition: push_controls.hpp:19
Simd< double > z
Definition: particles.hpp:19
constexpr long long int ps_was_inside
Definition: particle_stream.hpp:43
int ipc
Definition: push_controls.hpp:20
Definition: magnetic_field.F90:1
constexpr long long int ps_is_tracked
Definition: particle_stream.hpp:44