Degas2 physics Bug Fixes

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;

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<Omega_h::Real> 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,16 +426,14 @@ void read_input_parameters(int argc, char *const *argv,
 }
 
 void sample_initial_particle_energy(Kokkos::View<double *> 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
     // will give you a velocity vector (vx,vy,vz), a unit direction vector
     // (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<double *> 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);

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<o::Vector<3>, 4> &tet_verts) {
   return basis;
 }
 
+/*
 OMEGA_H_DEVICE o::Vector<3>
 barycentric2real(const o::Few<o::Vector<3>, 4> &tet_verts,
                  const o::Vector<4> &bary) {
@@ -1105,6 +1115,18 @@ barycentric2real(const o::Few<o::Vector<3>, 4> &tet_verts,
 
   return real_coords;
 }
+*/
+
+OMEGA_H_DEVICE o::Vector<3>
+barycentric2real(const o::Few<o::Vector<3>, 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<Omega_h::Real> 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<Kokkos::DefaultExecutionSpace> random_pool;
+  Kokkos::Random_XorShift64_Pool<Kokkos::DefaultExecutionSpace> 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);

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,22 +70,22 @@ 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) {
 		varl += (output(i).position[0] - l)*(output(i).position[0] - l);
 	}
 	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));
   }