update_particle_flags.hpp

Go to the documentation of this file.

#ifndef UPDATE_PARTICLE_FLAGS_HPP
#define UPDATE_PARTICLE_FLAGS_HPP

#include "particle_stream.hpp"
#include "magnetic_field.hpp"
#include "push_controls.hpp"

#ifdef ESC_PTL
// In the following subroutine, following flags need to be set.
//    - is_escaped 
//    - is_divertor  --> after sheath
//    - is_outboard  
//    - is_to_write : when is_escpaed or is_divertor set
//    - was_inside : set after it is used to check is_escaped
//    - do not consider the situation that is_escpaed and is_divertor happen together.
template<class Device>
KOKKOS_INLINE_FUNCTION void update_particle_flag1(const PushControls &push_controls, const MagneticField<Device> &magnetic_field, SimdParticles &part_tmp, Simd<double> &psi)
{
    Simd<long long int> itmp;
    Simd<bool> bt0, bt1, bt2, bt3; // temp. boolean

    for (int i_simd = 0; i_simd < SIMD_SIZE; i_simd++){
        // zero out above 25 bits. (8 bits for flags and 17 bits for esc_step)            
        part_tmp.flag[i_simd] = part_tmp.flag[i_simd] & ps_flag_step_mask;

        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
        //bt1[i_simd] = psi[i_simd] < XPOINT_PSI && x.z[i_simd] > XPOINT_Z ; // new position is inside
        bt1[i_simd] = (magnetic_field.equil.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)
        bt2[i_simd] = (!bt1[i_simd]) && flag_check(part_tmp.flag[i_simd], ps_was_inside);  // is_outside && was_inside --> escaped
        bt3[i_simd] = part_tmp.ph.r[i_simd] > magnetic_field.equil.axis_r ; // is_outboard
                    
        // set is_escpaed, is_outboard, is_to_write  (only when bt0 false)
        itmp[i_simd] = part_tmp.flag[i_simd];
        flag_set(itmp[i_simd], bt2[i_simd], ps_is_escaped); // set is_escpaed as bt2
        flag_set(itmp[i_simd], bt2[i_simd], ps_is_to_write); // set is_to_write (always when is_escaped true.)
        flag_set(itmp[i_simd], bt3[i_simd], ps_is_outboard); // set is_outboard as bt3
        flag_set(itmp[i_simd], bt1[i_simd], ps_was_inside); // set was_insde to future poistion

        // Particles contributing to diffusion coefficient analysis are marked with ps_is_tracked.
        // Only particles in the scrape-off layer contribute to this calculation.
        // First, particles crossing the separatrix into the scrape-off layer need to have the ps_was_traced flag set
        // Condition: ps_is_escaped == true (set above)
        flag_set(itmp[i_simd], bt2[i_simd] || flag_check(itmp[i_simd],ps_is_tracked), ps_is_tracked);  
        
    }
        
    //Set escaped time step -- separate loop due to if stetement. -- SIMD might not work.
    for (int i_simd = 0; i_simd < SIMD_SIZE; i_simd++){
        if(bt2[i_simd]) {
            //escaped flag set true
            long long int l_step = push_controls.gstep; // gstep in long long int data type.
            itmp[i_simd] = itmp[i_simd] & ps_flag_mask; // zero out above 8 bits. Clearing old value. 
            itmp[i_simd] +=  (l_step << ps_flag_bits);  // makes 8 zeros and adding to itmp.
            // if l_step is 131072 (2^17) or above, overflow happens. Currently XGC output supports up to 99,999 time step in output name.
        }
    }

    // update particle flag only when is_to_write == false
    for (int i_simd = 0; i_simd < SIMD_SIZE; i_simd++){
        part_tmp.flag[i_simd] = (bt0[i_simd]) ? (part_tmp.flag[i_simd]) : (itmp[i_simd]); 
    }

    return;
}

// For Gordon Bell 2024, the condition for tracked particles is changed.
// When a particle is inside of a flux shell (finite depth like 0.8 < psi < 0.81) at a certain time step, 
// is_to_write and is_tracked will be true. Otherwise, it will be false.
// This update will not be called for a several time steps.
// It only needs is_to_write
template<class Device>
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)
{
    Simd<long long int> itmp;
    Simd<bool> bt0; // temp. boolean

    if(push_controls.ipc==2 && push_controls.gstep%push_controls.particle_stream.reset_period==0){
        for (int i_simd = 0; i_simd < SIMD_SIZE; i_simd++){
            // zero out above 25 bits. (8 bits for flags and 17 bits for esc_step)            
            part_tmp.flag[i_simd] = part_tmp.flag[i_simd] & ps_flag_step_mask;


            bt0[i_simd] = (push_controls.particle_stream.pin < psi[i_simd]) && ( psi[i_simd]<push_controls.particle_stream.pout); 

            // set is_to_write
            itmp[i_simd] = part_tmp.flag[i_simd];
            flag_set(itmp[i_simd], bt0[i_simd], ps_is_to_write); // set is_to_write 
            flag_set(itmp[i_simd], bt0[i_simd], ps_is_tracked); // set is_to_write 

            // update particle flag 
            part_tmp.flag[i_simd] = itmp[i_simd];

        }
    }

    return;
}

KOKKOS_INLINE_FUNCTION void update_particle_flag2(Simd<bool> &not_in_triangle, Simd<double> &dw, SimdParticles &part_tmp)
{
    Simd<long long int> itmp;
    Simd<bool> bt0; // temp. boolean

    // update divertor flag (+ is_to_write) when the particle hit the divertor
    // only send the particle with ipc==2 loop.
    // can this be better to move after sheath_calculaton??
    for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
        
        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

        //beware that the following is called only when at_least_one is true
        bool dw_nonzero = not_in_triangle[i_simd] & (dw[i_simd]!=0.0);

        itmp[i_simd]=part_tmp.flag[i_simd];
        flag_set(itmp[i_simd], dw_nonzero, ps_is_divertor); // set is_divertor as dw_nonzero
        flag_set(itmp[i_simd], dw_nonzero, ps_is_to_write); // set is_to_write the same as is_to_write

        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
    }

    // Pack dw value to flag
    //beware that the following is called only when at_least_one is true
    // flag has zero bits above flag_mask (3 bytes + 1bit) at this stage.
    // dw can use 38 bits. dw * 2^nb1 will be saved as integer value. 
    // Let dw be between min_dw and max_dw
    // This should not exceed    2^(63-nb1-nb2)  
    constexpr double min_dw=-256.;
    constexpr double max_dw= 256.;

    // following factor1 and factor2 should be consistent unpacking fortran routine in 
    constexpr int nb1=ps_dw_factor_bits;  // nb1 is mutiplication factor to make floating point to integer
    constexpr int nb2=ps_flag_step_bits;  // nb2 is mutiplication factor to clear last 8 + 17 bits.
    constexpr int factor1 = 1u << nb1; // 2^28 , not more than 31
    constexpr int factor2=  1u << nb2; // 2^25

    // loop over all particles to utilize SIMD
    // dw should be zero for the particles that did not hit the divertor
    // zero is added at the end.
    // even though it is not zero, adding dw should be harmless.
    for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
        dw[i_simd] = min(max_dw, max(min_dw, dw[i_simd]));
        itmp[i_simd] = dw[i_simd]   * static_cast<double>(factor1); // itmp is the integer to place at first 38 bits.
        itmp[i_simd] = itmp[i_simd] * factor2; // multiplication is used to conserve sign, instead of '<< nb2'

        // 2022-03-03: Removed, RK
        // 2022-03-21: Revived.
        part_tmp.flag[i_simd]+=itmp[i_simd]; // add dw and flag. itmp should have zeros in the last 8 bit.
    }
    return;
}
#endif
#endif