Merge pull request #1450 from paulromano/fix-ace-isotropic

Data fix for isotropic elastic scattering / parallelize source site generation

include/openmc/timer.h

@@ -21,6 +21,7 @@ extern Timer time_finalize;
 extern Timer time_inactive;
 extern Timer time_initialize;
 extern Timer time_read_xs;
+extern Timer time_sample_source;
 extern Timer time_tallies;
 extern Timer time_total;
 extern Timer time_transport;

openmc/data/reaction.py

@@ -1073,11 +1073,15 @@ class Reaction(EqualityMixin):
         if i_reaction < ace.nxs[5] + 1:
             # Check if angular distribution data exist
             loc = int(ace.xss[ace.jxs[8] + i_reaction])
-            if loc <= 0:
-                # Angular distribution is either given as part of a product
-                # angle-energy distribution or is not given at all (in which
-                # case isotropic scattering is assumed)
+            if loc < 0:
+                # Angular distribution is given as part of a product
+                # angle-energy distribution
                 angle_dist = None
+            elif loc == 0:
+                # Angular distribution is isotropic
+                energy = [0.0, grid[-1]]
+                mu = Uniform(-1., 1.)
+                angle_dist = AngleDistribution(energy, [mu, mu])
             else:
                 angle_dist = AngleDistribution.from_ace(ace, ace.jxs[9], loc)
 

src/output.cpp

@@ -474,6 +474,8 @@ void print_runtime()
     show_time("Time synchronizing fission bank", time_bank.elapsed(), 1);
     show_time("Sampling source sites", time_bank_sample.elapsed(), 2);
     show_time("SEND/RECV source sites", time_bank_sendrecv.elapsed(), 2);
+  } else {
+    show_time("Time sampling source", time_sample_source.elapsed(), 1);
   }
   show_time("Time accumulating tallies", time_tallies.elapsed(), 1);
   show_time("Total time for finalization", time_finalize.elapsed());

src/physics.cpp

@@ -696,12 +696,12 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT,
   // Find speed of neutron in CM
   vel = v_n.norm();
 
-  // Sample scattering angle, checking if it is an ncorrelated angle-energy
-  // distribution
+  // Sample scattering angle, checking if angle distribution is present (assume
+  // isotropic otherwise)
   double mu_cm;
   auto& d = rx.products_[0].distribution_[0];
   auto d_ = dynamic_cast<UncorrelatedAngleEnergy*>(d.get());
-  if (d_) {
+  if (!d_->angle().empty()) {
     mu_cm = d_->angle().sample(p->E_, p->current_seed());
   } else {
     mu_cm = 2.0*prn(p->current_seed()) - 1.0;

src/simulation.cpp

@@ -462,7 +462,9 @@ void finalize_generation()
 
   } else if (settings::run_mode == RUN_MODE_FIXEDSOURCE) {
     // For fixed-source mode, we need to sample the external source
+    simulation::time_sample_source.start();
     fill_source_bank_fixedsource();
+    simulation::time_sample_source.stop();
   }
 }
 

src/source.cpp

@@ -328,6 +328,7 @@ void free_memory_source()
 void fill_source_bank_fixedsource()
 {
   if (settings::path_source.empty()) {
+    #pragma omp parallel for
     for (int64_t i = 0; i < simulation::work_per_rank; ++i) {
       // initialize random number seed
       int64_t id = (simulation::total_gen + overall_generation()) *

src/state_point.cpp

@@ -283,6 +283,8 @@ openmc_statepoint_write(const char* filename, bool* write_source)
       write_dataset(runtime_group, "synchronizing fission bank", time_bank.elapsed());
       write_dataset(runtime_group, "sampling source sites", time_bank_sample.elapsed());
       write_dataset(runtime_group, "SEND-RECV source sites", time_bank_sendrecv.elapsed());
+    } else {
+      write_dataset(runtime_group, "sampling source sites", time_sample_source.elapsed());
     }
     write_dataset(runtime_group, "accumulating tallies", time_tallies.elapsed());
     write_dataset(runtime_group, "total", time_total.elapsed());

src/timer.cpp

@@ -16,6 +16,7 @@ Timer time_finalize;
 Timer time_inactive;
 Timer time_initialize;
 Timer time_read_xs;
+Timer time_sample_source;
 Timer time_tallies;
 Timer time_total;
 Timer time_transport;
@@ -68,6 +69,7 @@ void reset_timers()
   simulation::time_inactive.reset();
   simulation::time_initialize.reset();
   simulation::time_read_xs.reset();
+  simulation::time_sample_source.reset();
   simulation::time_tallies.reset();
   simulation::time_total.reset();
   simulation::time_transport.reset();