From 5bffda07d5d8b08a18202e5e7fa5d8004b099a67 Mon Sep 17 00:00:00 2001
From: Ryan Thomas Wollaeger <wollaeger@lanl.gov>
Date: Wed, 29 Jan 2025 12:28:36 -0700
Subject: [PATCH] Do not divide particle energy weight by frequency.

+ Remove division by frequency in source particle creation.
+ Update singularity-opac commit to correct mean opacity accesses.
+ Have regressions fail if mean errors are NaN.
---
 external/singularity-opac | 2 +-
 src/jaybenne/sourcing.cpp | 5 +++--
 tst/regression_test.py    | 4 ++--
 3 files changed, 6 insertions(+), 5 deletions(-)

diff --git a/external/singularity-opac b/external/singularity-opac
index b61be52..b05c6e1 160000
--- a/external/singularity-opac
+++ b/external/singularity-opac
@@ -1 +1 @@
-Subproject commit b61be5226fee505d05ffeb6419408740792b51c7
+Subproject commit b05c6e12f003edae4cea38f6aa279c395aa11098
diff --git a/src/jaybenne/sourcing.cpp b/src/jaybenne/sourcing.cpp
index 3921bf5..733c4f3 100644
--- a/src/jaybenne/sourcing.cpp
+++ b/src/jaybenne/sourcing.cpp
@@ -223,6 +223,9 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) {
           ppack_i(b, ph::ijk(1), n) = j;
           ppack_i(b, ph::ijk(2), n) = k;
 
+          // Set energy weight
+          ppack_r(b, ph::weight(), n) = vmesh(b, fj::source_ew_per_cell(), k, j, i);
+
           // Sample position uniformly in space over cell
           // TODO(BRR) only valid for Cartesian
           ppack_r(b, swarm_position::x(), n) = xi + dx_i * (rng_gen.drand() - 0.5);
@@ -257,8 +260,6 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) {
             const Real nu = numind * std::exp((n + 0.5) * dlnu);
             ppack_r(b, ph::energy(), n) = hd * nu;
           }
-          ppack_r(b, ph::weight(), n) =
-              vmesh(b, fj::source_ew_per_cell(), k, j, i) / ppack_r(b, ph::energy(), n);
 
           if constexpr (ST == SourceType::emission) {
             dejbn -= ppack_r(b, ph::weight(), n);
diff --git a/tst/regression_test.py b/tst/regression_test.py
index 47aae00..981757d 100644
--- a/tst/regression_test.py
+++ b/tst/regression_test.py
@@ -412,10 +412,10 @@ def analytic_comparison(
     print(f"Max fractional error:           {max_frac_error:.2e}")
 
     if args.comparison == "mean":
-        if mean_frac_error > tolerance:
+        if mean_frac_error > tolerance or np.isnan(mean_frac_error):
             success = False
     elif args.comparison == "weighted_mean":
-        if mean_frac_error_weighted > tolerance:
+        if mean_frac_error_weighted > tolerance or np.isnan(mean_frac_error):
             success = False
 
     if args.visualize: