gistfile1.txt

// in tracer_diag.cpp
    plasma.for_all_nonadiabatic_species([&](Species<DeviceType>& species){
        // This is the syntax to access the Cabana AoSoA
        // See example code in https://github.com/ECP-copa/Cabana/blob/master/core/unit_test/tstAoSoA.hpp
        auto phase = Cabana::slice<0>(species.particles);
        auto tracer = Cabana::slice<2>(species.particles);
        auto gid = Cabana::slice<3>(species.particles);
        auto flag = Cabana::slice<4>(species.particles);
        
        written_ptls = 0;
        log_msg.str("");
        log_msg << "\n tracer_diag - calc loop. Particles = " << species.n_ptl << "\n\n";
        checkpoint(log_msg.str(), Chkpt::NoBarrier);

        // Zero out flag counters
        ctr_is_to_write = 0;
        ctr_is_in_init = 0;
        ctr_is_written = 0;
        ctr_is_escaped = 0;
        ctr_is_divertor = 0;
        ctr_is_outboard = 0;
        ctr_was_inside = 0;
        ctr_is_traced = 0;

        // Species.n_ptl is the number of particles per MPI rank.
        for (int p = 0; p < species.n_ptl; p++) { // Loop over particles in the current species
            if( (flag(p) & tracer_is_to_write) > 0) {
                ctr_is_to_write++;
            }
            if( (flag(p) & tracer_is_in_init) > 0) {
                ctr_is_in_init++;
            }
            if( (flag(p) & tracer_is_written) > 0) {
                ctr_is_written++;
            }
            if( (flag(p) & tracer_is_escaped) > 0) {
                ctr_is_escaped++;
            }
            if( (flag(p) & tracer_is_divertor) > 0) {
                ctr_is_divertor++;
            }
            if( (flag(p) & tracer_is_outboard) > 0) {
                ctr_is_outboard++;
            }
            if( (flag(p) & tracer_was_inside) > 0) {
                ctr_was_inside++;
            }
            if( (flag(p) & tracer_is_traced) > 0) {
                ctr_is_traced++;  
            }

            if(gid(p) == 98445)
            {
                std::cout << "gid = " << gid(p) << ", R = " << phase(p, 0) << ", Z = " << phase(p, 1).z << ", flag = " << flag(p) << std::endl;
            }



////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// boundaries.tpp


    if(sml.ipc==2 ){
        // part_tmp.flag is the same with part_one.flag here -- data copied before.
        // In the following loop, 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.
        for (int i_simd = 0; i_simd < SIMD_SIZE; i_simd++){
            if(part_tmp.gid[i_simd] == 98445){
                   printf("\n\n Before logic block: particle gid %llu: R=%f Z=%f psi=%f is_in_region1=%d bt0=%d bt1=%d bt2=%d bt3=%d itmp=%d flag=%d\n\n", 
                            part_tmp.gid[i_simd], x.r[i_simd], x.z[i_simd], psi[i_simd],
                            magnetic_field.equil.is_in_region_1(x.r[i_simd],x.z[i_simd],psi[i_simd]),
                            bt0[i_simd], bt1[i_simd], bt2[i_simd], bt3[i_simd], itmp[i_simd], part_tmp.flag[i_simd]);
            }
            // zero out abvoe 25 bits. (8 bits for flags and 17 bits for esc_step)            
            part_tmp.flag[i_simd] = part_tmp.flag[i_simd] & tracer_flag_step_mask;

            /**** This is the original code */
            bt0[i_simd] = (part_tmp.flag[i_simd] & tracer_is_to_write) > 0; // it needs to be written (determined previously), do not change any flag
            bt1[i_simd] = !(magnetic_field.equil.is_in_region_1(x.r[i_simd],x.z[i_simd],psi[i_simd])); //  new position is not region 1 ( = outside)
            bt2[i_simd] = bt1[i_simd] && ((part_tmp.flag[i_simd] & tracer_was_inside) > 0);  // is_outside && was_inside --> escaped
            bt3[i_simd] = x.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];
            itmp[i_simd] = (bt2[i_simd])? ( itmp[i_simd] | tracer_is_escaped ) : ( itmp[i_simd] & ~tracer_is_escaped  ); // set is_escpaed
            itmp[i_simd] = (bt2[i_simd])? ( itmp[i_simd] | tracer_is_to_write) : ( itmp[i_simd] & ~tracer_is_to_write ); // set is_to_write
            itmp[i_simd] = (bt3[i_simd])? ( itmp[i_simd] | tracer_is_outboard) : ( itmp[i_simd] & ~tracer_is_outboard ); // set is_outboard
            itmp[i_simd] =(!bt1[i_simd])? ( itmp[i_simd] | tracer_was_inside ) : ( itmp[i_simd] & ~tracer_was_inside  ); // update was_inside to future position
            

            part_tmp.flag[i_simd] = (bt0[i_simd]) ? (part_tmp.flag[i_simd]) : (itmp[i_simd]); // update it only is_to_write == false
            if(part_tmp.gid[i_simd] == 98445){
                   printf("\n\n After logic block: particle gid %llu: R=%f Z=%f psi=%f is_in_region1=%d bt0=%d bt1=%d bt2=%d bt3=%d itmp=%d flag=%d\n\n", 
                            part_tmp.gid[i_simd], x.r[i_simd], x.z[i_simd], psi[i_simd],
                            magnetic_field.equil.is_in_region_1(x.r[i_simd],x.z[i_simd],psi[i_simd]),
                            bt0[i_simd], bt1[i_simd], bt2[i_simd], bt3[i_simd], itmp[i_simd], part_tmp.flag[i_simd]);
               }
}