Created
March 2, 2022 22:19
-
-
Save rkube/32303bb58dbdc8d584ecc985c5c83522 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
// 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]); | |
} | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment