From baee04637411bf351f796168f013363fc6e6afa9 Mon Sep 17 00:00:00 2001 From: gpaynter5004 Date: Thu, 24 Jul 2025 16:10:12 -0400 Subject: [PATCH 01/11] Update degas2-main.cpp Fix bulk_flow_velocity size --- src/physics/degas2/degas2-main.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/physics/degas2/degas2-main.cpp b/src/physics/degas2/degas2-main.cpp index eaf45d3..76f61b8 100644 --- a/src/physics/degas2/degas2-main.cpp +++ b/src/physics/degas2/degas2-main.cpp @@ -335,7 +335,7 @@ void get_field_values(Omega_h::Reals centroids, Fields &fields) { fields.ion_temperature.resize(num_elements, 1.0); fields.electron_density.resize(num_elements, 1.0); fields.ion_density.resize(num_elements, 1.0); - fields.bulk_flow_velocity.resize(num_elements, 0.0); + fields.bulk_flow_velocity.resize(3*num_elements, 0.0); // TODO: Add your code here to retrieve the actual field values // Write centroid coords to CSV to be read by python script From c7fc5e193bc928bc6ecfaa57b551f3b25c95d117 Mon Sep 17 00:00:00 2001 From: Gabe Date: Fri, 25 Jul 2025 12:10:04 -0400 Subject: [PATCH 02/11] Set the alpha to actually update --- src/physics/degas2/DG2Physics.h | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/physics/degas2/DG2Physics.h b/src/physics/degas2/DG2Physics.h index 31be246..71374d7 100644 --- a/src/physics/degas2/DG2Physics.h +++ b/src/physics/degas2/DG2Physics.h @@ -155,6 +155,9 @@ class DG2Physics { particle_info.position[0] += l*particle_info.direction[0]; particle_info.position[1] += l*particle_info.direction[1]; particle_info.position[2] += l*particle_info.direction[2]; + + double mp {938.27e6/(3e10*3e10)}; //eV/c^2 = eV*s^2/cm^2 + particle_info.alpha = Kokkos::sqrt(mp/(2.0*particle_energy(particle_info.particle_index))); } // collision event @@ -231,6 +234,9 @@ class DG2Physics { } } particle_info.weight = new_weight; + + //This definition may be redundant since it is also in the next location function + particle_info.alpha = Kokkos::sqrt(mp/(2.0*particle_energy(particle_info.particle_index))); random_pool.free_state(rand_gen); From b5852dc6250278049cd94f3c91244a9be3e27015 Mon Sep 17 00:00:00 2001 From: Gabe Date: Fri, 25 Jul 2025 13:18:56 -0400 Subject: [PATCH 03/11] Fix formatting --- src/physics/degas2/DG2Physics.h | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/physics/degas2/DG2Physics.h b/src/physics/degas2/DG2Physics.h index 71374d7..733d2a8 100644 --- a/src/physics/degas2/DG2Physics.h +++ b/src/physics/degas2/DG2Physics.h @@ -238,9 +238,7 @@ class DG2Physics { //This definition may be redundant since it is also in the next location function particle_info.alpha = Kokkos::sqrt(mp/(2.0*particle_energy(particle_info.particle_index))); random_pool.free_state(rand_gen); - - -} + } //To here random_pool_t random_pool; DG2CrossSection cross_section; From 952647dbf79738694f190588adc5dbd5ff99f223 Mon Sep 17 00:00:00 2001 From: "Gabe J. Paynter" Date: Fri, 25 Jul 2025 18:16:48 -0400 Subject: [PATCH 04/11] Bug fix: random_pool, uniform source distribution, NEEDS TO FIX: illegal memory access --- src/physics/degas2/degas2-main.cpp | 12 +++-- src/pumitallyopenmc/pumitally_impl.cpp | 62 +++++++++++++++++++------- 2 files changed, 55 insertions(+), 19 deletions(-) diff --git a/src/physics/degas2/degas2-main.cpp b/src/physics/degas2/degas2-main.cpp index 76f61b8..63581db 100644 --- a/src/physics/degas2/degas2-main.cpp +++ b/src/physics/degas2/degas2-main.cpp @@ -72,7 +72,7 @@ int main(int argc, char *argv[]) { // Initialize PUMI-Tally and Read Fields auto pumi_tally = pumiinopenmc::PumiTallyImpl( - input_params.mesh_name, input_params.num_particles, argc, argv); + input_params.mesh_name, input_params.num_particles, argc, argv, input_params.source_distribution); auto &mesh = pumi_tally.full_mesh_; Omega_h::Write centroids(mesh.nelems() * 3); get_centroids(mesh, centroids); @@ -135,11 +135,14 @@ void transport(pumiinopenmc::PumiTallyImpl &pumi_tally, DG2Physics &physics, auto last_exit = pumi_tally.p_pumi_particle_at_elem_boundary_handler->last_exit_; + auto alpha = pumi_tally.p_pumi_particle_at_elem_boundary_handler->alpha_; auto get_new_destination = PS_LAMBDA(const int &e, const int &pid, const int &mask) { + if (mask > 0) { // FIXME: check if the particle is at destination or at // the boundary + // ParticleInfo particle_info; particle_info.position[0] = particle_dest(pid, 0); particle_info.position[1] = particle_dest(pid, 1); @@ -177,6 +180,7 @@ void transport(pumiinopenmc::PumiTallyImpl &pumi_tally, DG2Physics &physics, particle_group(pid) = particle_info.energy_group; alpha[pid] = particle_info.alpha; + } }; pumipic::parallel_for(pumi_tally.pumipic_ptcls.get(), get_new_destination, @@ -411,8 +415,7 @@ void read_input_parameters(int argc, char *const *argv, } void sample_initial_particle_energy(Kokkos::View energy_array) { - random_pool_t randomPool; - + random_pool_t randomPool(0); auto sample_energy = OMEGA_H_LAMBDA(int i) { // Basically, feed this 4 uniformly generated random numbers, x1, x2, y1, // y2, on the interval (0,1) and it @@ -420,7 +423,6 @@ void sample_initial_particle_energy(Kokkos::View energy_array) { // (directionx, directiony, directionz), the alpha value used in the tally, // and the energy of the particle that was sampled. All from a gas at // temperature temp. - double mp{938.27e6 / (3e10 * 3e10)}; // eV/c^2 = eV*s^2/cm^2. Necessary // constant for the distribution @@ -465,6 +467,8 @@ void sample_initial_particle_energy(Kokkos::View energy_array) { // computations. double particle_energy = 0.5 * mp * mag_v * mag_v; energy_array(i) = particle_energy; + + }; Kokkos::parallel_for("sample_initial_particle_energy", energy_array.size(), sample_energy); diff --git a/src/pumitallyopenmc/pumitally_impl.cpp b/src/pumitallyopenmc/pumitally_impl.cpp index 9aae074..d49a3c4 100644 --- a/src/pumitallyopenmc/pumitally_impl.cpp +++ b/src/pumitallyopenmc/pumitally_impl.cpp @@ -592,18 +592,31 @@ void PumiParticleAtElemBoundary::updatePrevXPoint( void PumiParticleAtElemBoundary::updatePrevXPoint(PPPS *ptcls) { // todo add checks of size auto prev_xpoints_l = prev_xpoint_; + printf("\n Size of previous xpoints %d\n", prev_xpoints_l.size()); + printf("\nNumber of particles: %d\n", ptcls->nPtcls()); OMEGA_H_CHECK_PRINTF( ptcls->capacity() * 3 == prev_xpoints_l.size(), "Error: prev_xpoints_s are not size properly capacity %d size %d\n", ptcls->capacity(), prev_xpoints_l.size()); + printf("\n Omega_h Check Ran \n"); auto xpoints = ptcls->get<0>(); + printf("\n** xpoints defined **\n"); auto update = PS_LAMBDA(const auto &e, const auto &pid, const auto &mask) { + printf("\n** Inside lambda ** \n"); + if (mask>0 && pid<10) { + + printf("\n ***Copying xpoints***\n"); prev_xpoints_l[pid * 3 + 0] = xpoints(pid, 0); prev_xpoints_l[pid * 3 + 1] = xpoints(pid, 1); prev_xpoints_l[pid * 3 + 2] = xpoints(pid, 2); + } }; + printf("\n ** After the loop ** \n"); pumipic::parallel_for(ptcls, update, "update previous xpoints from origin points"); + printf("\nBefore fence \n"); + Kokkos::fence(); + printf("\n after fence \n"); } void PumiParticleAtElemBoundary::evaluateFlux( @@ -618,6 +631,7 @@ void PumiParticleAtElemBoundary::evaluateFlux( auto p_wgt = ptcls->get<4>(); auto p_groups = ptcls->get<5>(); auto xpoints_l = xpoints; // todo shouldn't need it, so remove + auto alpha_l = alpha_; auto evaluate_flux = PS_LAMBDA(const int &e, const int &pid, const int &mask) { @@ -648,7 +662,7 @@ void PumiParticleAtElemBoundary::evaluateFlux( orig[0], orig[1], orig[2], dest[0], dest[1], dest[2]); } - Omega_h::Real contribution = segment_length * p_wgt(pid) * alpha_[pid]; + Omega_h::Real contribution = segment_length * p_wgt(pid) * alpha_l[pid]; int group = p_groups(pid); OMEGA_H_CHECK_PRINTF( @@ -873,6 +887,7 @@ PumiPIC particle structure size to %d\n", n_particles); void start_pumi_particles_in_0th_element(Omega_h::Mesh &mesh, pumiinopenmc::PPPS *ptcls) { // find the centroid of the 0th element + printf("\n***Started 0th Element***\n"); const auto &coords = mesh.coords(); const auto &tet2node = mesh.ask_down(Omega_h::REGION, Omega_h::VERT).ab2b; @@ -943,12 +958,16 @@ void PumiTallyImpl::create_and_initialize_pumi_particle_structure( init_loc(pid, 0) = device_pos_buffer_l[pid * 3 + 0]; init_loc(pid, 1) = device_pos_buffer_l[pid * 3 + 1]; init_loc(pid, 2) = device_pos_buffer_l[pid * 3 + 2]; + if (pid < 10){ + printf("\n PID %d in element %d: (%f, %f, %f)\n",pid, e, init_loc(pid,0), init_loc(pid,1), init_loc(pid,2)); + } } }; pumipic::parallel_for(pumipic_ptcls.get(), copy_initial_positions, "copy initial positions from device buffer"); } if (source_dist == SourceDistribution::ZERO) { + printf("\**nSource Distribution 0**\n"); start_pumi_particles_in_0th_element(*mesh, pumipic_ptcls.get()); } p_pumi_particle_at_elem_boundary_handler = @@ -1090,6 +1109,7 @@ barycentric_basis(const o::Few, 4> &tet_verts) { return basis; } +/* OMEGA_H_DEVICE o::Vector<3> barycentric2real(const o::Few, 4> &tet_verts, const o::Vector<4> &bary) { @@ -1105,6 +1125,18 @@ barycentric2real(const o::Few, 4> &tet_verts, return real_coords; } +*/ + +OMEGA_H_DEVICE o::Vector<3> +barycentric2real(const o::Few, 4>& tet_verts, +const o::Vector<4>& bary) { +o::Vector<3> real_coords {0,0,0}; +for (int i = 0; i < 4; ++i) { +real_coords += bary[i] * tet_verts[i]; +} +return real_coords; +} + void initialize_uniform_source(Omega_h::Mesh &mesh, Omega_h::Write particle_positions, @@ -1135,7 +1167,7 @@ void initialize_uniform_source(Omega_h::Mesh &mesh, const auto cells2nodes = mesh.ask_down(o::REGION, o::VERT).ab2b; const auto coords = mesh.coords(); - Kokkos::Random_XorShift64_Pool random_pool; + Kokkos::Random_XorShift64_Pool random_pool(0); auto set_initial_positions = OMEGA_H_LAMBDA(const int &e) { auto pid_start = cumulative_particles[e]; auto pid_end = cumulative_particles[e + 1]; @@ -1144,19 +1176,19 @@ void initialize_uniform_source(Omega_h::Mesh &mesh, for (Omega_h::LO pid = pid_start; pid < pid_end; ++pid) { auto gen = random_pool.get_state(); - o::Vector<4> random_bcc{0.0, 0.0, 0.0, 0.0}; - random_bcc[0] = gen.drand(0.0, 1.0); - random_bcc[1] = gen.drand(0.0, 1.0); - random_bcc[2] = gen.drand(0.0, 1.0); - o::Real complimentary0 = 1.0 - random_bcc[0]; - o::Real complimentary1 = 1.0 - random_bcc[1]; - o::Real complimentary2 = 1.0 - random_bcc[2]; - - bool more_than_one = random_bcc[0] + random_bcc[1] + random_bcc[2] > 1.0; - random_bcc[0] = more_than_one ? complimentary0 : random_bcc[0]; - random_bcc[1] = more_than_one ? complimentary1 : random_bcc[1]; - random_bcc[2] = more_than_one ? complimentary2 : random_bcc[2]; - random_bcc[2] = 1.0 - random_bcc[0] - random_bcc[1] - random_bcc[2]; + o::Real r1 = gen.drand(0.0, 1.0); + o::Real r2 = gen.drand(0.0, 1.0); + o::Real r3 = gen.drand(0.0, 1.0); + + r1 = Kokkos::pow(r1, 1.0 / 3.0); + r2 = Kokkos::sqrt(r2); + o::Real a = 1.0 - r1; + o::Real b = r1 * (1.0 - r2); + o::Real c = r1 * r2 * (1.0 - r3); + o::Real d = r1 * r2 * r3; + + o::Vector<4> random_bcc{a, b, c, d}; + random_pool.free_state(gen); auto verts = o::gather_verts<4>(cells2nodes, e); From 6c09ac0d6e53f890eb97ee8bb6fb98ea1ef2f2a2 Mon Sep 17 00:00:00 2001 From: Gabe Paynter Date: Tue, 29 Jul 2025 20:32:58 -0400 Subject: [PATCH 05/11] Fix illegal memory access error. NEED TO FIX: particles deleting --- src/physics/degas2/DG2Physics.h | 12 ++++++------ src/physics/degas2/degas2-main.cpp | 18 +++++++++++------- src/pumitallyopenmc/pumitally_impl.cpp | 20 +++++--------------- 3 files changed, 22 insertions(+), 28 deletions(-) diff --git a/src/physics/degas2/DG2Physics.h b/src/physics/degas2/DG2Physics.h index 71374d7..29d5116 100644 --- a/src/physics/degas2/DG2Physics.h +++ b/src/physics/degas2/DG2Physics.h @@ -151,7 +151,7 @@ class DG2Physics { double sigma_cx = charge_exchange_cross_section(particle_info, field_info); //Generate distance and move particle - double l =-Kokkos::log(x)/(field_info.electron_density*sigma_ion+field_info.ion_density*sigma_cx); //cm. n in cm^-3 + double l =-Kokkos::log(x)/(field_info.electron_density*sigma_ion+field_info.ion_density*sigma_cx)*0.01; //m. n in cm^-3 particle_info.position[0] += l*particle_info.direction[0]; particle_info.position[1] += l*particle_info.direction[1]; particle_info.position[2] += l*particle_info.direction[2]; @@ -176,12 +176,12 @@ class DG2Physics { //Compute 3 Maxwellian (Gaussian) distributed velocities (cm/s) double mp {938.27e6/(3e10*3e10)}; //eV/c^2 = eV*s^2/cm^2 double max_sigma_cx = 3.8e-14; - + bool rejection_test = false; auto old_mag_v = Kokkos::sqrt(2*particle_energy(particle_info.particle_index)/mp); - double vx; - double vy; - double vz; + double vx {}; + double vy {}; + double vz {}; //Loops this until it passes the rejection test while (!rejection_test) { //Generate random numbers for the velocity @@ -213,7 +213,7 @@ class DG2Physics { particle_info.direction[0] = vx/mag_v; particle_info.direction[1] = vy/mag_v; particle_info.direction[2] = vz/mag_v; - + particle_energy(particle_info.particle_index) = 0.5*mp*mag_v*mag_v; //particle_info.energy_group = particle_energy(particle_info.particle_index); //Temporary for debugging //Adjust Weights diff --git a/src/physics/degas2/degas2-main.cpp b/src/physics/degas2/degas2-main.cpp index 63581db..2def0be 100644 --- a/src/physics/degas2/degas2-main.cpp +++ b/src/physics/degas2/degas2-main.cpp @@ -166,12 +166,17 @@ void transport(pumiinopenmc::PumiTallyImpl &pumi_tally, DG2Physics &physics, field_info.bulk_flow_velocity[0] = bulk_flow_velocity[e * 3 + 0]; field_info.bulk_flow_velocity[1] = bulk_flow_velocity[e * 3 + 1]; field_info.bulk_flow_velocity[2] = bulk_flow_velocity[e * 3 + 2]; - - if (last_exit[pid] == -1) { // reached destination - physics.collide_particle(particle_info, field_info); - } + + if (iter != 0) { //Ensure last_exit[pid] doesn't run the first time + if (last_exit[pid] == -1) { // reached destination + physics.collide_particle(particle_info, field_info); + } + } + else { + physics.collide_particle(particle_info, field_info); + } physics.sample_collision_distance(particle_info, field_info); - + // Update particle position and direction particle_dest(pid, 0) = particle_info.position[0]; particle_dest(pid, 1) = particle_info.position[1]; @@ -180,12 +185,11 @@ void transport(pumiinopenmc::PumiTallyImpl &pumi_tally, DG2Physics &physics, particle_group(pid) = particle_info.energy_group; alpha[pid] = particle_info.alpha; - + printf("\n#### PID: %d, iter: %d, at position (%f,%f,%f) ####\n", pid, iter, particle_info.position[0], particle_info.position[1], particle_info.position[2]); } }; pumipic::parallel_for(pumi_tally.pumipic_ptcls.get(), get_new_destination, "get new destination"); - pumi_tally.search_and_rebuild( false, true); // for now, always rebuild the pp structure Kokkos::fence(); diff --git a/src/pumitallyopenmc/pumitally_impl.cpp b/src/pumitallyopenmc/pumitally_impl.cpp index d49a3c4..42430d0 100644 --- a/src/pumitallyopenmc/pumitally_impl.cpp +++ b/src/pumitallyopenmc/pumitally_impl.cpp @@ -592,31 +592,21 @@ void PumiParticleAtElemBoundary::updatePrevXPoint( void PumiParticleAtElemBoundary::updatePrevXPoint(PPPS *ptcls) { // todo add checks of size auto prev_xpoints_l = prev_xpoint_; - printf("\n Size of previous xpoints %d\n", prev_xpoints_l.size()); - printf("\nNumber of particles: %d\n", ptcls->nPtcls()); OMEGA_H_CHECK_PRINTF( ptcls->capacity() * 3 == prev_xpoints_l.size(), "Error: prev_xpoints_s are not size properly capacity %d size %d\n", ptcls->capacity(), prev_xpoints_l.size()); - printf("\n Omega_h Check Ran \n"); auto xpoints = ptcls->get<0>(); - printf("\n** xpoints defined **\n"); auto update = PS_LAMBDA(const auto &e, const auto &pid, const auto &mask) { - printf("\n** Inside lambda ** \n"); if (mask>0 && pid<10) { - printf("\n ***Copying xpoints***\n"); prev_xpoints_l[pid * 3 + 0] = xpoints(pid, 0); prev_xpoints_l[pid * 3 + 1] = xpoints(pid, 1); prev_xpoints_l[pid * 3 + 2] = xpoints(pid, 2); } }; - printf("\n ** After the loop ** \n"); pumipic::parallel_for(ptcls, update, - "update previous xpoints from origin points"); - printf("\nBefore fence \n"); - Kokkos::fence(); - printf("\n after fence \n"); + "update previous xpoints from origin points"); } void PumiParticleAtElemBoundary::evaluateFlux( @@ -958,16 +948,16 @@ void PumiTallyImpl::create_and_initialize_pumi_particle_structure( init_loc(pid, 0) = device_pos_buffer_l[pid * 3 + 0]; init_loc(pid, 1) = device_pos_buffer_l[pid * 3 + 1]; init_loc(pid, 2) = device_pos_buffer_l[pid * 3 + 2]; - if (pid < 10){ - printf("\n PID %d in element %d: (%f, %f, %f)\n",pid, e, init_loc(pid,0), init_loc(pid,1), init_loc(pid,2)); - } +// if (pid < 10){ +// printf("\n PID %d in element %d: (%f, %f, %f)\n",pid, e, init_loc(pid,0), init_loc(pid,1), init_loc(pid,2)); +// } } }; pumipic::parallel_for(pumipic_ptcls.get(), copy_initial_positions, "copy initial positions from device buffer"); } if (source_dist == SourceDistribution::ZERO) { - printf("\**nSource Distribution 0**\n"); + printf("\n**Source Distribution 0**\n"); start_pumi_particles_in_0th_element(*mesh, pumipic_ptcls.get()); } p_pumi_particle_at_elem_boundary_handler = From f13992a890dc0d902471a4f6ca3b4bc7f74829ec Mon Sep 17 00:00:00 2001 From: Gabe Paynter Date: Fri, 1 Aug 2025 14:29:08 -0400 Subject: [PATCH 06/11] Fix particle deleting due to improper definitions of igin and destination position. NEED TO FIX: new issue of particles deleting --- src/physics/degas2/DG2Physics.h | 9 ++++++--- src/physics/degas2/degas2-main.cpp | 23 +++++++++++++++-------- 2 files changed, 21 insertions(+), 11 deletions(-) diff --git a/src/physics/degas2/DG2Physics.h b/src/physics/degas2/DG2Physics.h index a1c160c..2e5f2d9 100644 --- a/src/physics/degas2/DG2Physics.h +++ b/src/physics/degas2/DG2Physics.h @@ -151,13 +151,16 @@ class DG2Physics { double sigma_cx = charge_exchange_cross_section(particle_info, field_info); //Generate distance and move particle +// printf("\n$$$ Before moving particle [%d], position: (%f, %f, %f) $$$\n", particle_info.particle_index, particle_info.position[0], particle_info.position[1], +// particle_info.position[2]); double l =-Kokkos::log(x)/(field_info.electron_density*sigma_ion+field_info.ion_density*sigma_cx)*0.01; //m. n in cm^-3 particle_info.position[0] += l*particle_info.direction[0]; particle_info.position[1] += l*particle_info.direction[1]; particle_info.position[2] += l*particle_info.direction[2]; - +// printf("\n$$$ After moving particle [%d], position (%f, %f, %f). Direction was (%f, %f, %f) with l = %f $$$\n", particle_info.particle_index, particle_info.position[0], +// particle_info.position[1], particle_info.position[2], particle_info.direction[0], particle_info.direction[1], particle_info.direction[2], l); double mp {938.27e6/(3e10*3e10)}; //eV/c^2 = eV*s^2/cm^2 - particle_info.alpha = Kokkos::sqrt(mp/(2.0*particle_energy(particle_info.particle_index))); + particle_info.alpha = Kokkos::sqrt(mp/(2.0*particle_energy(particle_info.particle_index)))*100; //1/(m/s) } // collision event @@ -236,7 +239,7 @@ class DG2Physics { particle_info.weight = new_weight; //This definition may be redundant since it is also in the next location function - particle_info.alpha = Kokkos::sqrt(mp/(2.0*particle_energy(particle_info.particle_index))); + particle_info.alpha = Kokkos::sqrt(mp/(2.0*particle_energy(particle_info.particle_index)))*100; //1/(m/s) random_pool.free_state(rand_gen); } //To here diff --git a/src/physics/degas2/degas2-main.cpp b/src/physics/degas2/degas2-main.cpp index 2def0be..755b15e 100644 --- a/src/physics/degas2/degas2-main.cpp +++ b/src/physics/degas2/degas2-main.cpp @@ -136,21 +136,23 @@ void transport(pumiinopenmc::PumiTallyImpl &pumi_tally, DG2Physics &physics, auto last_exit = pumi_tally.p_pumi_particle_at_elem_boundary_handler->last_exit_; + printf("\n### iter = %d ###\n",iter); + auto alpha = pumi_tally.p_pumi_particle_at_elem_boundary_handler->alpha_; auto get_new_destination = PS_LAMBDA(const int &e, const int &pid, const int &mask) { - + if (mask > 0) { // FIXME: check if the particle is at destination or at // the boundary // ParticleInfo particle_info; - particle_info.position[0] = particle_dest(pid, 0); - particle_info.position[1] = particle_dest(pid, 1); - particle_info.position[2] = particle_dest(pid, 2); + particle_info.position[0] = particle_orig(pid, 0); + particle_info.position[1] = particle_orig(pid, 1); + particle_info.position[2] = particle_orig(pid, 2); auto direction = Omega_h::normalize(Omega_h::Vector<3>{ - particle_info.position[0] - particle_orig(pid, 0), - particle_info.position[1] - particle_orig(pid, 1), - particle_info.position[2] - particle_orig(pid, 2)}); + particle_info.position[0] - particle_dest(pid, 0), + particle_info.position[1] - particle_dest(pid, 1), + particle_info.position[2] - particle_dest(pid, 2)}); particle_info.direction[0] = direction[0]; particle_info.direction[1] = direction[1]; particle_info.direction[2] = direction[2]; @@ -167,6 +169,10 @@ void transport(pumiinopenmc::PumiTallyImpl &pumi_tally, DG2Physics &physics, field_info.bulk_flow_velocity[1] = bulk_flow_velocity[e * 3 + 1]; field_info.bulk_flow_velocity[2] = bulk_flow_velocity[e * 3 + 2]; + +// printf("\n---- Direction (pid = %d) before collision: (%f, %f, %f). Original pos: (%f, %f, %f), pos (%f, %f, %f) ---- \n", pid, direction[0], direction[1], direction[2], +// particle_orig(pid,0), particle_orig(pid,1), particle_orig(pid,2), particle_info.position[0], particle_info.position[1], particle_info.position[2]); + if (iter != 0) { //Ensure last_exit[pid] doesn't run the first time if (last_exit[pid] == -1) { // reached destination physics.collide_particle(particle_info, field_info); @@ -185,7 +191,8 @@ void transport(pumiinopenmc::PumiTallyImpl &pumi_tally, DG2Physics &physics, particle_group(pid) = particle_info.energy_group; alpha[pid] = particle_info.alpha; - printf("\n#### PID: %d, iter: %d, at position (%f,%f,%f) ####\n", pid, iter, particle_info.position[0], particle_info.position[1], particle_info.position[2]); +// printf("\n#### PID: %d, iter: %d, at position (%f,%f,%f). Direction was: (%f, %f, %f) ####\n", pid, iter, particle_info.position[0], +// particle_info.position[1], particle_info.position[2], particle_info.direction[0], particle_info.direction[1], particle_info.direction[2]); } }; pumipic::parallel_for(pumi_tally.pumipic_ptcls.get(), get_new_destination, From e66f8a862622a6d25ac0abccfdb7a3a59247e055 Mon Sep 17 00:00:00 2001 From: gpaynter5004 Date: Fri, 8 Aug 2025 11:54:37 -0400 Subject: [PATCH 07/11] Re-enable the print for the test test_degas2_physics.cpp --- test/test_degas2_physics.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/test/test_degas2_physics.cpp b/test/test_degas2_physics.cpp index 28154e1..03d853d 100644 --- a/test/test_degas2_physics.cpp +++ b/test/test_degas2_physics.cpp @@ -61,7 +61,7 @@ TEST_CASE("Test Degas2 Physics Functions"){ }); auto output = create_mirror_view(particles); Kokkos::deep_copy(output, particles); - //std::ofstream outfile("Log.txt"); + std::ofstream outfile("Log.txt"); //These 4 functions should return the commented value for SEED=12345 @@ -70,7 +70,7 @@ TEST_CASE("Test Degas2 Physics Functions"){ l += output(i).position[0]; } l /= numParticles; - //outfile << "Average Distance (cm): " << l << std::endl; //2.41 + outfile << "Average Distance (cm): " << l << std::endl; //2.41 double varl {0}; for (int i=0; i < numParticles; ++i) { @@ -78,14 +78,14 @@ TEST_CASE("Test Degas2 Physics Functions"){ } varl /= (numParticles - 1); double sdl {sqrt(varl)}; - //outfile << "Standard Deviation of Distance (cm): " << sdl << std::endl; //2.38 + outfile << "Standard Deviation of Distance (cm): " << sdl << std::endl; //2.38 double ux {0}; for (int i=0; i < numParticles; ++i) { ux += output(i).direction[0]; } ux /= numParticles; - //outfile << "Mean x Direction: " << ux << std::endl; //0.0159 + outfile << "Mean x Direction: " << ux << std::endl; //0.0159 double varux; for (int i=0; i < numParticles; ++i) { @@ -94,7 +94,7 @@ TEST_CASE("Test Degas2 Physics Functions"){ varux /= (numParticles - 1); double sdux {sqrt(varux)}; - /* + outfile << "Standard Deviation of Mean x Direction: " << sdux << std::endl; //0.579 @@ -102,7 +102,7 @@ TEST_CASE("Test Degas2 Physics Functions"){ for (int i = 0; i < numParticles; ++i) { outfile << i << ","; - outfile << output(i).energy_group << ","; + outfile << (1.0/output(i).alpha) << ","; outfile << output(i).weight << ","; outfile << output(i).position[0] << ","; outfile << output(i).position[1] << ","; @@ -111,7 +111,7 @@ TEST_CASE("Test Degas2 Physics Functions"){ outfile << output(i).direction[1] << ","; outfile << output(i).direction[2] << std::endl; } - */ + REQUIRE_THAT(l, Catch::Matchers::WithinAbs(2.40,.04)); REQUIRE_THAT(sdl, Catch::Matchers::WithinAbs(2.40,.04)); REQUIRE_THAT(ux, Catch::Matchers::WithinAbs(0.0,.03)); From 2ad7a739874a54abf2abde69d9bec168df5be358 Mon Sep 17 00:00:00 2001 From: gpaynter5004 Date: Fri, 8 Aug 2025 11:57:58 -0400 Subject: [PATCH 08/11] Correct rejection sampling DG2Physics.h --- src/physics/degas2/DG2Physics.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/physics/degas2/DG2Physics.h b/src/physics/degas2/DG2Physics.h index 2e5f2d9..58b8d05 100644 --- a/src/physics/degas2/DG2Physics.h +++ b/src/physics/degas2/DG2Physics.h @@ -178,7 +178,7 @@ class DG2Physics { //Compute New Direction and Energy and set particle info //Compute 3 Maxwellian (Gaussian) distributed velocities (cm/s) double mp {938.27e6/(3e10*3e10)}; //eV/c^2 = eV*s^2/cm^2 - double max_sigma_cx = 3.8e-14; + double max_rate_cx = 1.4e-7; bool rejection_test = false; auto old_mag_v = Kokkos::sqrt(2*particle_energy(particle_info.particle_index)/mp); @@ -206,7 +206,7 @@ class DG2Physics { auto mag_v2 = rel_vx*rel_vx + rel_vy*rel_vy + rel_vz*rel_vz; //Generate random number and compare to sigma/sigma_max - if (rand_gen.drand(0., 1.) < analytic_cross_section_CE(0.5*mp*mag_v2)/max_sigma_cx) { + if (rand_gen.drand(0., 1.) < mag_v2*analytic_cross_section_CE(0.5*mp*mag_v2)/max_rate_cx) { rejection_test = true; } } From 29fa14efc7f615be3be700db65b56516f68ebb2d Mon Sep 17 00:00:00 2001 From: gpaynter5004 Date: Fri, 8 Aug 2025 12:17:59 -0400 Subject: [PATCH 09/11] Correct units in test case test_degas2_physics.cpp --- test/test_degas2_physics.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_degas2_physics.cpp b/test/test_degas2_physics.cpp index 03d853d..9b107b2 100644 --- a/test/test_degas2_physics.cpp +++ b/test/test_degas2_physics.cpp @@ -112,8 +112,8 @@ TEST_CASE("Test Degas2 Physics Functions"){ outfile << output(i).direction[2] << std::endl; } - REQUIRE_THAT(l, Catch::Matchers::WithinAbs(2.40,.04)); - REQUIRE_THAT(sdl, Catch::Matchers::WithinAbs(2.40,.04)); + REQUIRE_THAT(l, Catch::Matchers::WithinAbs(0.0240,0.0004)); + REQUIRE_THAT(sdl, Catch::Matchers::WithinAbs(0.0240,0.0004)); REQUIRE_THAT(ux, Catch::Matchers::WithinAbs(0.0,.03)); REQUIRE_THAT(sdux, Catch::Matchers::WithinAbs(0.577,.03)); } From 7c601900bb8f83649e386c91a3f8d196e2c6ebbe Mon Sep 17 00:00:00 2001 From: gpaynter5004 Date: Fri, 8 Aug 2025 12:21:04 -0400 Subject: [PATCH 10/11] Correct units in log test_degas2_physics.cpp --- test/test_degas2_physics.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_degas2_physics.cpp b/test/test_degas2_physics.cpp index 9b107b2..daef4d7 100644 --- a/test/test_degas2_physics.cpp +++ b/test/test_degas2_physics.cpp @@ -98,7 +98,7 @@ TEST_CASE("Test Degas2 Physics Functions"){ outfile << "Standard Deviation of Mean x Direction: " << sdux << std::endl; //0.579 - outfile << "Particle #,Energy(eV),Weight,X(cm),Y(cm),Z(cm),X_Dir,Y_Dir,Z_Dir" << std::endl; + outfile << "Particle #,Velocity(m/s),Weight,X(cm),Y(cm),Z(cm),X_Dir,Y_Dir,Z_Dir" << std::endl; for (int i = 0; i < numParticles; ++i) { outfile << i << ","; From b401232bf274730549ae78dbb326cb87768ee407 Mon Sep 17 00:00:00 2001 From: gpaynter5004 Date: Fri, 8 Aug 2025 12:46:53 -0400 Subject: [PATCH 11/11] Bug fix rejection sampling DG2Physics.h --- src/physics/degas2/DG2Physics.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/physics/degas2/DG2Physics.h b/src/physics/degas2/DG2Physics.h index 58b8d05..e1549ae 100644 --- a/src/physics/degas2/DG2Physics.h +++ b/src/physics/degas2/DG2Physics.h @@ -206,7 +206,7 @@ class DG2Physics { auto mag_v2 = rel_vx*rel_vx + rel_vy*rel_vy + rel_vz*rel_vz; //Generate random number and compare to sigma/sigma_max - if (rand_gen.drand(0., 1.) < mag_v2*analytic_cross_section_CE(0.5*mp*mag_v2)/max_rate_cx) { + if (rand_gen.drand(0., 1.) < Kokkos::sqrt(mag_v2)*analytic_cross_section_CE(0.5*mp*mag_v2)/max_rate_cx) { rejection_test = true; } }