Merge branch 'develop' into IFP_beta_nuc

Author: GuySten

.github/workflows/ci.yml

@@ -1,4 +1,4 @@
-name: CI
+name: Tests and Coverage
 
 on:
   # allows us to run workflows manually
@@ -21,7 +21,25 @@ env:
   GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
 
 jobs:
+  filter-changes:
+    runs-on: ubuntu-latest
+    outputs:
+      source_changed: ${{ steps.filter.outputs.source_changed }}
+    steps:
+      - name: Check out the repository
+        uses: actions/checkout@v4
+      - name: Examine changed files
+        id: filter
+        uses: dorny/paths-filter@668c092af3649c4b664c54e4b704aa46782f6f7c # latest master commit, not released yet
+        with:
+          filters: |
+            source_changed:
+              - '!docs/**'
+              - '!**/*.md'
+          predicate-quantifier: 'every'
   main:
+    needs: filter-changes
+    if: ${{ needs.filter-changes.outputs.source_changed == 'true' }}
     runs-on: ubuntu-22.04
     strategy:
       matrix:
@@ -205,12 +223,33 @@ jobs:
           flag-name: C++ and Python
           path-to-lcov: coverage.lcov
 
-  finish:
-    needs: main
+  coverage:
+    needs: [filter-changes, main]
+    if: ${{ always() }}
     runs-on: ubuntu-latest
     steps:
       - name: Coveralls Finished
+        if: ${{ needs.filter-changes.outputs.source_changed == 'true' }}
         uses: coverallsapp/github-action@v2
         with:
           github-token: ${{ secrets.GITHUB_TOKEN }}
           parallel-finished: true
+
+  ci-pass:
+    needs: [filter-changes, main, coverage]
+    name: Check CI status
+    if: ${{ always() }}
+    runs-on: ubuntu-latest
+    steps:
+      - name: Check CI status
+        run: |
+          if [[ "${{ needs.filter-changes.outputs.source_changed }}" == "false" ]]; then
+            echo "Documentation-only change - CI skipped successfully"
+            exit 0
+          fi
+          if [[ "${{ needs.main.result }}" == "success" && "${{ needs.coverage.result }}" == "success" ]]; then
+            echo "CI passed"
+            exit 0
+          fi
+          echo "CI failed"
+          exit 1

AGENTS.md

@@ -31,6 +31,24 @@ C++17/Python codebase where:
 - **Depletion**: `openmc.deplete` implements burnup via operator-splitting with various integrators (Predictor, CECM, etc.)
 - **Nuclear Data**: `openmc.data` provides programmatic access to nuclear data files (ENDF, ACE, HDF5)
 
+## Git Branching Workflow
+
+OpenMC uses a git flow branching model with two primary branches:
+
+- **`develop` branch**: The main development branch where all ongoing development takes place. This is the **primary branch against which pull requests are submitted and merged**. This branch is not guaranteed to be stable and may contain work-in-progress features.
+- **`master` branch**: The stable release branch containing the latest stable release of OpenMC. This branch only receives merges from `develop` when the development team decides a release should occur.
+
+### Instructions for Code Review
+
+When analyzing code changes on a feature or bugfix branch (e.g., when a user asks "what do you think of these changes?"), **compare the branch changes against `develop`, not `master`**. Pull requests are submitted to merge into `develop`, so differences relative to `develop` represent the actual proposed changes. Comparing against `master` will include unrelated changes from other features that have already been merged to `develop`.
+
+### Workflow for contributors
+
+1. Create a feature/bugfix branch off `develop`
+2. Make changes and commit to the feature branch
+3. Open a pull request to merge the feature branch into `develop`
+4. A committer reviews and merges the PR into `develop`
+
 ## Critical Build & Test Workflows
 
 ### Build Dependencies

include/openmc/material.h

@@ -14,6 +14,7 @@
 #include "openmc/memory.h" // for unique_ptr
 #include "openmc/ncrystal_interface.h"
 #include "openmc/particle.h"
+#include "openmc/settings.h"
 #include "openmc/vector.h"
 
 namespace openmc {
@@ -110,9 +111,12 @@ class Material {
   //! \return Density in [atom/b-cm]
   double density() const { return density_; }
 
-  //! Get density in [g/cm^3]
+  //! Get density in [g/cm^3].
   //! \return Density in [g/cm^3]
-  double density_gpcc() const { return density_gpcc_; }
+  double density_gpcc() const
+  {
+    return settings::run_CE ? density_gpcc_ : density();
+  }
 
   //! Get charge density in [e/b-cm]
   //! \return Charge density in [e/b-cm]

include/openmc/random_ray/source_region.h

@@ -146,6 +146,7 @@ class SourceRegionHandle {
 
   // Scalar fields
   int* material_;
+  double* density_mult_;
   int* is_small_;
   int* n_hits_;
   int* birthday_;
@@ -195,6 +196,9 @@ class SourceRegionHandle {
   int& material() { return *material_; }
   const int material() const { return *material_; }
 
+  double& density_mult() { return *density_mult_; }
+  const double density_mult() const { return *density_mult_; }
+
   int& is_small() { return *is_small_; }
   const int is_small() const { return *is_small_; }
 
@@ -316,7 +320,9 @@ class SourceRegion {
   //---------------------------------------
   // Scalar fields
 
-  int material_ {0}; //!< Index in openmc::model::materials array
+  int material_ {0};          //!< Index in openmc::model::materials array
+  double density_mult_ {1.0}; //!< A density multiplier queried from the cell
+                              //!< corresponding to the source region.
   OpenMPMutex lock_;
   double volume_ {
     0.0}; //!< Volume (computed from the sum of ray crossing lengths)
@@ -394,6 +400,9 @@ class SourceRegionContainer {
   int& material(int64_t sr) { return material_[sr]; }
   const int material(int64_t sr) const { return material_[sr]; }
 
+  double& density_mult(int64_t sr) { return density_mult_[sr]; }
+  const double density_mult(int64_t sr) const { return density_mult_[sr]; }
+
   int& is_small(int64_t sr) { return is_small_[sr]; }
   const int is_small(int64_t sr) const { return is_small_[sr]; }
 
@@ -625,6 +634,7 @@ class SourceRegionContainer {
 
   // SoA storage for scalar fields (one item per source region)
   vector<int> material_;
+  vector<double> density_mult_;
   vector<int> is_small_;
   vector<int> n_hits_;
   vector<int> mesh_;

src/random_ray/flat_source_domain.cpp

@@ -109,18 +109,21 @@ void FlatSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
 
   // Add scattering + fission source
   int material = srh.material();
+  double density_mult = srh.density_mult();
   if (material != MATERIAL_VOID) {
     double inverse_k_eff = 1.0 / k_eff_;
     for (int g_out = 0; g_out < negroups_; g_out++) {
-      double sigma_t = sigma_t_[material * negroups_ + g_out];
+      double sigma_t = sigma_t_[material * negroups_ + g_out] * density_mult;
       double scatter_source = 0.0;
       double fission_source = 0.0;
 
       for (int g_in = 0; g_in < negroups_; g_in++) {
         double scalar_flux = srh.scalar_flux_old(g_in);
-        double sigma_s =
-          sigma_s_[material * negroups_ * negroups_ + g_out * negroups_ + g_in];
-        double nu_sigma_f = nu_sigma_f_[material * negroups_ + g_in];
+        double sigma_s = sigma_s_[material * negroups_ * negroups_ +
+                                  g_out * negroups_ + g_in] *
+                         density_mult;
+        double nu_sigma_f =
+          nu_sigma_f_[material * negroups_ + g_in] * density_mult;
         double chi = chi_[material * negroups_ + g_out];
 
         scatter_source += sigma_s * scalar_flux;
@@ -198,7 +201,8 @@ void FlatSourceDomain::set_flux_to_flux_plus_source(
         source_regions_.volume_sq(sr);
     }
   } else {
-    double sigma_t = sigma_t_[source_regions_.material(sr) * negroups_ + g];
+    double sigma_t = sigma_t_[source_regions_.material(sr) * negroups_ + g] *
+                     source_regions_.density_mult(sr);
     source_regions_.scalar_flux_new(sr, g) /= (sigma_t * volume);
     source_regions_.scalar_flux_new(sr, g) += source_regions_.source(sr, g);
   }
@@ -332,7 +336,8 @@ void FlatSourceDomain::compute_k_eff()
     double sr_fission_source_new = 0;
 
     for (int g = 0; g < negroups_; g++) {
-      double nu_sigma_f = nu_sigma_f_[material * negroups_ + g];
+      double nu_sigma_f = nu_sigma_f_[material * negroups_ + g] *
+                          source_regions_.density_mult(sr);
       sr_fission_source_old +=
         nu_sigma_f * source_regions_.scalar_flux_old(sr, g);
       sr_fission_source_new +=
@@ -562,7 +567,8 @@ double FlatSourceDomain::compute_fixed_source_normalization_factor() const
       // to get the total source strength in the expected units.
       double sigma_t = 1.0;
       if (material != MATERIAL_VOID) {
-        sigma_t = sigma_t_[material * negroups_ + g];
+        sigma_t =
+          sigma_t_[material * negroups_ + g] * source_regions_.density_mult(sr);
       }
       simulation_external_source_strength +=
         source_regions_.external_source(sr, g) * sigma_t * volume;
@@ -618,7 +624,9 @@ void FlatSourceDomain::random_ray_tally()
     // source strength.
     double volume = source_regions_.volume(sr) * simulation_volume_;
 
-    double material = source_regions_.material(sr);
+    int material = source_regions_.material(sr);
+    double density_mult = source_regions_.density_mult(sr);
+
     for (int g = 0; g < negroups_; g++) {
       double flux =
         source_regions_.scalar_flux_new(sr, g) * source_normalization_factor;
@@ -634,19 +642,22 @@ void FlatSourceDomain::random_ray_tally()
 
         case SCORE_TOTAL:
           if (material != MATERIAL_VOID) {
-            score = flux * volume * sigma_t_[material * negroups_ + g];
+            score =
+              flux * volume * sigma_t_[material * negroups_ + g] * density_mult;
           }
           break;
 
         case SCORE_FISSION:
           if (material != MATERIAL_VOID) {
-            score = flux * volume * sigma_f_[material * negroups_ + g];
+            score =
+              flux * volume * sigma_f_[material * negroups_ + g] * density_mult;
           }
           break;
 
         case SCORE_NU_FISSION:
           if (material != MATERIAL_VOID) {
-            score = flux * volume * nu_sigma_f_[material * negroups_ + g];
+            score = flux * volume * nu_sigma_f_[material * negroups_ + g] *
+                    density_mult;
           }
           break;
 
@@ -913,7 +924,8 @@ void FlatSourceDomain::output_to_vtk() const
             for (int g = 0; g < negroups_; g++) {
               int64_t source_element = fsr * negroups_ + g;
               float flux = evaluate_flux_at_point(voxel_positions[i], fsr, g);
-              double sigma_f = sigma_f_[mat * negroups_ + g];
+              double sigma_f = sigma_f_[mat * negroups_ + g] *
+                               source_regions_.density_mult(fsr);
               total_fission += sigma_f * flux;
             }
           }
@@ -934,7 +946,8 @@ void FlatSourceDomain::output_to_vtk() const
             // multiply it back to get the true external source.
             double sigma_t = 1.0;
             if (mat != MATERIAL_VOID) {
-              sigma_t = sigma_t_[mat * negroups_ + g];
+              sigma_t = sigma_t_[mat * negroups_ + g] *
+                        source_regions_.density_mult(fsr);
             }
             total_external += source_regions_.external_source(fsr, g) * sigma_t;
           }
@@ -1244,7 +1257,8 @@ void FlatSourceDomain::set_adjoint_sources()
       continue;
     }
     for (int g = 0; g < negroups_; g++) {
-      double sigma_t = sigma_t_[material * negroups_ + g];
+      double sigma_t =
+        sigma_t_[material * negroups_ + g] * source_regions_.density_mult(sr);
       source_regions_.external_source(sr, g) /= sigma_t;
     }
   }
@@ -1495,6 +1509,8 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
 
   handle.material() = material;
 
+  handle.density_mult() = cell.density_mult(gs.cell_instance());
+
   // Store the mesh index (if any) assigned to this source region
   handle.mesh() = mesh_idx;
 
@@ -1523,7 +1539,8 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
     // Divide external source term by sigma_t
     if (material != C_NONE) {
       for (int g = 0; g < negroups_; g++) {
-        double sigma_t = sigma_t_[material * negroups_ + g];
+        double sigma_t =
+          sigma_t_[material * negroups_ + g] * handle.density_mult();
         handle.external_source(g) /= sigma_t;
       }
     }
@@ -1598,16 +1615,18 @@ void FlatSourceDomain::apply_transport_stabilization()
 #pragma omp parallel for
   for (int64_t sr = 0; sr < n_source_regions(); sr++) {
     int material = source_regions_.material(sr);
+    double density_mult = source_regions_.density_mult(sr);
     if (material == MATERIAL_VOID) {
       continue;
     }
     for (int g = 0; g < negroups_; g++) {
       // Only apply stabilization if the diagonal (in-group) scattering XS is
       // negative
       double sigma_s =
-        sigma_s_[material * negroups_ * negroups_ + g * negroups_ + g];
+        sigma_s_[material * negroups_ * negroups_ + g * negroups_ + g] *
+        density_mult;
       if (sigma_s < 0.0) {
-        double sigma_t = sigma_t_[material * negroups_ + g];
+        double sigma_t = sigma_t_[material * negroups_ + g] * density_mult;
         double phi_new = source_regions_.scalar_flux_new(sr, g);
         double phi_old = source_regions_.scalar_flux_old(sr, g);
 

src/random_ray/linear_source_domain.cpp

@@ -43,12 +43,13 @@ void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
 
   // Add scattering + fission source
   int material = srh.material();
+  double density_mult = srh.density_mult();
   if (material != MATERIAL_VOID) {
     double inverse_k_eff = 1.0 / k_eff_;
     MomentMatrix invM = srh.mom_matrix().inverse();
 
     for (int g_out = 0; g_out < negroups_; g_out++) {
-      double sigma_t = sigma_t_[material * negroups_ + g_out];
+      double sigma_t = sigma_t_[material * negroups_ + g_out] * density_mult;
 
       double scatter_flat = 0.0f;
       double fission_flat = 0.0f;
@@ -61,9 +62,11 @@ void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
         MomentArray flux_linear = srh.flux_moments_old(g_in);
 
         // Handles for cross sections
-        double sigma_s =
-          sigma_s_[material * negroups_ * negroups_ + g_out * negroups_ + g_in];
-        double nu_sigma_f = nu_sigma_f_[material * negroups_ + g_in];
+        double sigma_s = sigma_s_[material * negroups_ * negroups_ +
+                                  g_out * negroups_ + g_in] *
+                         density_mult;
+        double nu_sigma_f =
+          nu_sigma_f_[material * negroups_ + g_in] * density_mult;
         double chi = chi_[material * negroups_ + g_out];
 
         // Compute source terms for flat and linear components of the flux

src/random_ray/random_ray.cpp

@@ -435,7 +435,8 @@ void RandomRay::attenuate_flux_flat_source(
 
   // MOC incoming flux attenuation + source contribution/attenuation equation
   for (int g = 0; g < negroups_; g++) {
-    float sigma_t = domain_->sigma_t_[material * negroups_ + g];
+    float sigma_t =
+      domain_->sigma_t_[material * negroups_ + g] * srh.density_mult();
     float tau = sigma_t * distance;
     float exponential = cjosey_exponential(tau); // exponential = 1 - exp(-tau)
     float new_delta_psi = (angular_flux_[g] - srh.source(g)) * exponential;
@@ -558,7 +559,8 @@ void RandomRay::attenuate_flux_linear_source(
   for (int g = 0; g < negroups_; g++) {
 
     // Compute tau, the optical thickness of the ray segment
-    float sigma_t = domain_->sigma_t_[material * negroups_ + g];
+    float sigma_t =
+      domain_->sigma_t_[material * negroups_ + g] * srh.density_mult();
     float tau = sigma_t * distance;
 
     // If tau is very small, set it to zero to avoid numerical issues.