store progeny energy in SourceSite E_last

Author: GuySten

include/openmc/particle_data.h

@@ -44,6 +44,7 @@ struct SourceSite {
   Position r;
   Direction u;
   double E;
+  double E_last;
   double time {0.0};
   double wgt {1.0};
   int delayed_group {0};

openmc/lib/core.py

@@ -22,6 +22,7 @@ class _SourceSite(Structure):
     _fields_ = [('r', c_double*3),
                 ('u', c_double*3),
                 ('E', c_double),
+                ('E_last', c_double),                
                 ('time', c_double),
                 ('wgt', c_double),
                 ('delayed_group', c_int),

src/initialize.cpp

@@ -157,26 +157,27 @@ void initialize_mpi(MPI_Comm intracomm)
 
   // Create bank datatype
   SourceSite b;
-  MPI_Aint disp[11];
+  MPI_Aint disp[12];
   MPI_Get_address(&b.r, &disp[0]);
   MPI_Get_address(&b.u, &disp[1]);
   MPI_Get_address(&b.E, &disp[2]);
-  MPI_Get_address(&b.time, &disp[3]);
-  MPI_Get_address(&b.wgt, &disp[4]);
-  MPI_Get_address(&b.delayed_group, &disp[5]);
-  MPI_Get_address(&b.surf_id, &disp[6]);
-  MPI_Get_address(&b.particle, &disp[7]);
-  MPI_Get_address(&b.parent_nuclide, &disp[8]);
-  MPI_Get_address(&b.parent_id, &disp[9]);
-  MPI_Get_address(&b.progeny_id, &disp[10]);
-  for (int i = 10; i >= 0; --i) {
+  MPI_Get_address(&b.E_last, &disp[3]);
+  MPI_Get_address(&b.time, &disp[4]);
+  MPI_Get_address(&b.wgt, &disp[5]);
+  MPI_Get_address(&b.delayed_group, &disp[6]);
+  MPI_Get_address(&b.surf_id, &disp[7]);
+  MPI_Get_address(&b.particle, &disp[8]);
+  MPI_Get_address(&b.parent_nuclide, &disp[9]);
+  MPI_Get_address(&b.parent_id, &disp[10]);
+  MPI_Get_address(&b.progeny_id, &disp[11]);
+  for (int i = 11; i >= 0; --i) {
     disp[i] -= disp[0];
   }
 
-  int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1};
+  int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
   MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE,
     MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG};
-  MPI_Type_create_struct(11, blocks, disp, types, &mpi::source_site);
+  MPI_Type_create_struct(12, blocks, disp, types, &mpi::source_site);
   MPI_Type_commit(&mpi::source_site);
 }
 #endif // OPENMC_MPI

src/mesh.cpp

@@ -321,6 +321,7 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
 
     SourceSite site;
     site.E = 1.0;
+    site.E_last = 1.0;
     site.particle = ParticleType::neutron;
 
     for (int axis = 0; axis < 3; ++axis) {

src/particle.cpp

@@ -132,6 +132,7 @@ void Particle::from_source(const SourceSite* src, ParticleType particle)
   r_last_current() = src->r;
   r_last() = src->r;
   u_last() = src->u;
+  E_last() = src->E_last();
   if (settings::run_CE) {
     E() = src->E;
     g() = 0;
@@ -140,7 +141,6 @@ void Particle::from_source(const SourceSite* src, ParticleType particle)
     g_last() = static_cast<int>(src->E);
     E() = data::mg.energy_bin_avg_[g()];
   }
-  E_last() = E();
   type_last() = particle;
   time() = src->time;
   time_last() = src->time;

src/physics.cpp

@@ -209,6 +209,7 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
     site.r = p.r();
     site.particle = ParticleType::neutron;
     site.time = p.time();
+    site.E_last = p.E();
     site.wgt = 1. / weight;
     site.parent_id = p.id();
     site.progeny_id = p.n_progeny()++;

src/physics_mg.cpp

@@ -138,6 +138,7 @@ void create_fission_sites(Particle& p)
     site.r = p.r();
     site.particle = ParticleType::neutron;
     site.time = p.time();
+    site.E_last = p.E();
     site.wgt = 1. / weight;
     site.parent_id = p.id();
     site.progeny_id = p.n_progeny()++;

src/random_ray/random_ray.cpp

@@ -796,6 +796,7 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain)
   }
 
   site.E = 0.0;
+  site.E_last = 0.0;
   this->from_source(&site, site.particle);
 
   // Locate ray