diff --git a/src/physics/degas2/DG2Physics.h b/src/physics/degas2/DG2Physics.h index 31be246..e1549ae 100644 --- a/src/physics/degas2/DG2Physics.h +++ b/src/physics/degas2/DG2Physics.h @@ -151,10 +151,16 @@ 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 +// 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)))*100; //1/(m/s) } // collision event @@ -172,13 +178,13 @@ 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); - 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 @@ -200,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.) < Kokkos::sqrt(mag_v2)*analytic_cross_section_CE(0.5*mp*mag_v2)/max_rate_cx) { rejection_test = true; } } @@ -210,7 +216,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 @@ -231,10 +237,11 @@ class DG2Physics { } } particle_info.weight = new_weight; - random_pool.free_state(rand_gen); - -} + //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)))*100; //1/(m/s) + random_pool.free_state(rand_gen); + } //To here random_pool_t random_pool; DG2CrossSection cross_section; diff --git a/src/physics/degas2/degas2-main.cpp b/src/physics/degas2/degas2-main.cpp index eaf45d3..755b15e 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,19 +135,24 @@ 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]; @@ -163,12 +168,21 @@ 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); - } + + +// 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); + } + } + 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]; @@ -177,11 +191,12 @@ 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). 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, "get new destination"); - pumi_tally.search_and_rebuild( false, true); // for now, always rebuild the pp structure Kokkos::fence(); @@ -335,7 +350,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 @@ -411,8 +426,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 +434,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 +478,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..42430d0 100644 --- a/src/pumitallyopenmc/pumitally_impl.cpp +++ b/src/pumitallyopenmc/pumitally_impl.cpp @@ -598,12 +598,15 @@ void PumiParticleAtElemBoundary::updatePrevXPoint(PPPS *ptcls) { ptcls->capacity(), prev_xpoints_l.size()); auto xpoints = ptcls->get<0>(); auto update = PS_LAMBDA(const auto &e, const auto &pid, const auto &mask) { + if (mask>0 && pid<10) { + 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); + } }; pumipic::parallel_for(ptcls, update, - "update previous xpoints from origin points"); + "update previous xpoints from origin points"); } void PumiParticleAtElemBoundary::evaluateFlux( @@ -618,6 +621,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 +652,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 +877,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 +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)); +// } } }; pumipic::parallel_for(pumipic_ptcls.get(), copy_initial_positions, "copy initial positions from device buffer"); } if (source_dist == SourceDistribution::ZERO) { + printf("\n**Source Distribution 0**\n"); start_pumi_particles_in_0th_element(*mesh, pumipic_ptcls.get()); } p_pumi_particle_at_elem_boundary_handler = @@ -1090,6 +1099,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 +1115,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 +1157,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 +1166,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); diff --git a/test/test_degas2_physics.cpp b/test/test_degas2_physics.cpp index 28154e1..daef4d7 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,15 +94,15 @@ 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 - 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 << ","; - 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,9 +111,9 @@ 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(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)); }