lanl · RyanWollaeger · Jul 30, 2026 · Mar 19, 2026 · Apr 7, 2026 · Apr 8, 2026
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -141,3 +141,5 @@ jobs:
           ./stepdiff_smr.py --executable $JAYBENNE_EXECUTABLE \
                             --input ../inputs/stepdiff_smr_hybrid.in --use_mpiexec \
                             --mpi_nthreads 8 --mpi_oversubscribe
+          ./stepdiff_mg_plaw.py --executable $JAYBENNE_EXECUTABLE \
+                            --input ../inputs/stepdiff_mg_plaw_ddmc.in --use_mpiexec
diff --git a/analysis/plot.py b/analysis/plot.py
@@ -23,7 +23,16 @@
 
 # Plot a 1D profile of a variable
 def plot_1d(
-    fig, ax, filenames, variable_name, draw_meshblocks, vmin, vmax, coords, scale
+    fig,
+    ax,
+    filenames,
+    variable_name,
+    draw_meshblocks,
+    vmin,
+    vmax,
+    data_vbnd,
+    coords,
+    scale,
 ):
 
     for filename in filenames:
@@ -39,6 +48,11 @@ def plot_1d(
         if scale == "log":
             variable = np.log10(variable)
 
+        if data_vbnd:
+            for b in range(dump.NumBlocks):
+                vmin = min(vmin, np.min(variable[b, idx_k, idx_j, :]))
+                vmax = max(vmax, np.max(variable[b, idx_k, idx_j, :]))
+
         for b in range(dump.NumBlocks):
             ax.plot(dump.xc[b, idx_k, idx_j, :], variable[b, idx_k, idx_j, :])
 
@@ -54,6 +68,7 @@ def plot_2d(
     draw_meshblocks,
     vmin,
     vmax,
+    data_vbnd,
     coords,
     scale,
 ):
@@ -70,6 +85,11 @@ def plot_2d(
     if scale == "log":
         variable = np.log10(variable)
 
+    if data_vbnd:
+        for b in range(dump.NumBlocks):
+            vmin = min(vmin, np.min(variable[b, idx_k, :, :]))
+            vmax = max(vmax, np.max(variable[b, idx_k, :, :]))
+
     for b in range(dump.NumBlocks):
         ax.pcolormesh(
             dump.xc[b, idx_k, :, :],
@@ -161,6 +181,9 @@ def plot_2d(
     parser.add_argument(
         "--vmax", type=float, default=0, help="Maximum value of colorbar"
     )
+    parser.add_argument(
+        "--data_vbnd", action="store_true", help="Use min/max data for colorbar"
+    )
     parser.add_argument(
         "--scale",
         type=str,
@@ -187,6 +210,7 @@ def plot_2d(
             args.meshblocks,
             args.vmin,
             args.vmax,
+            args.data_vbnd,
             args.coords,
             args.scale,
         )
@@ -201,6 +225,7 @@ def plot_2d(
             args.meshblocks,
             args.vmin,
             args.vmax,
+            args.data_vbnd,
             args.coords,
             args.scale,
         )

diff --git a/external/singularity-opac b/external/singularity-opac
diff --git a/inputs/inf.in b/inputs/inf.in
@@ -59,15 +59,19 @@ seed = 349857
 
 <mcblock>
 frequency_type = gray
-opacity_model = constant
-opacity_constant_value = 1.0 # cm^2/g
-scattering_model = constant
-scattering_constant_value = 1.0e5
 specific_heat = 1.0e8 # erg/K/g
 initial_density = 1.0 # g cm^-3
 initial_temperature = 1.0e5 # K
 initial_radiation = thermal
 
+<mcblock/absorption>
+opacity_model = constant
+constant_value = 1.0 # cm^2/g
+
+<mcblock/scattering>
+opacity_model = constant
+constant_value = 1.0e5
+
 <parthenon/output0>
 file_type = hdf5
 dt = 0.1

diff --git a/inputs/inf_mg.in b/inputs/inf_mg.in
@@ -54,14 +54,18 @@ seed = 349851
 
 <mcblock>
 frequency_type = multigroup
-opacity_model = ep_bremss
-scattering_model = constant
-scattering_constant_value = 1e-1
 specific_heat = 1.0e8 # erg/K/g
 initial_density = 1.0e-5 # g cm^-3
 initial_temperature = 1.0e5 # K
 initial_radiation = thermal
 
+<mcblock/absorption>
+opacity_model = ep_bremss
+
+<mcblock/scattering>
+scattering_model = constant
+scattering_constant_value = 1e-1
+
 <parthenon/output0>
 file_type = hdf5
 dt = 0.1

diff --git a/inputs/inf_stiff.in b/inputs/inf_stiff.in
@@ -54,15 +54,19 @@ seed = 349856
 
 <mcblock>
 frequency_type = gray
-opacity_model = constant
-opacity_constant_value = 1000.0 # cm^2/g
-scattering_model = constant
-scattering_constant_value = 0.0 #1.0e5
 specific_heat = 1.0e7 # erg/K/g
 initial_density = 1.0 # g cm^-3
 initial_temperature = 1.0e5 # K
 initial_radiation = thermal
 
+<mcblock/absorption>
+opacity_model = constant
+constant_value = 1000.0 # cm^2/g
+
+<mcblock/scattering>
+opacity_model = constant
+constant_value = 0.0 #1.0e5
+
 <parthenon/output0>
 file_type = hdf5
 dt = 1e-11

diff --git a/inputs/stepdiff.in b/inputs/stepdiff.in
@@ -66,14 +66,18 @@ seed = 349857
 
 <mcblock>
 frequency_type = gray
-opacity_model = none
-scattering_model = constant
-scattering_constant_value = 1.0e3
 specific_heat = 1.0e8 # erg/K/g
 initial_density = 1.0 # g cm^-3
 initial_temperature = 1.0e5 # K
 initial_radiation = thermal
 
+<mcblock/absorption>
+opacity_model = none
+
+<mcblock/scattering>
+opacity_model = constant
+constant_value = 1.0e3
+
 <parthenon/output0>
 file_type = hdf5
 dt = 3.335641e-10

diff --git a/inputs/stepdiff_ddmc.in b/inputs/stepdiff_ddmc.in
@@ -67,14 +67,18 @@ seed = 349857
 
 <mcblock>
 frequency_type = gray
-opacity_model = none
-scattering_model = constant
-scattering_constant_value = 1.0e3
 specific_heat = 1.0e8 # erg/K/g
 initial_density = 1.0 # g cm^-3
 initial_temperature = 1.0e5 # K
 initial_radiation = thermal
 
+<mcblock/absorption>
+opacity_model = none
+
+<mcblock/scattering>
+opacity_model = constant
+constant_value = 1.0e3
+
 <parthenon/output0>
 file_type = hdf5
 dt = 3.335641e-10

diff --git a/inputs/stepdiff_mg_plaw_ddmc.in b/inputs/stepdiff_mg_plaw_ddmc.in
@@ -0,0 +1,96 @@
+# ========================================================================================
+#  (C) (or copyright) 2023-2024. Triad National Security, LLC. All rights reserved.
+#
+#  This program was produced under U.S. Government contract 89233218CNA000001 for Los
+#  Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+#  for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+#  in the program are reserved by Triad National Security, LLC, and the U.S. Department
+#  of Energy/National Nuclear Security Administration. The Government is granted for
+#  itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+#  license in this material to reproduce, prepare derivative works, distribute copies to
+#  the public, perform publicly and display publicly, and to permit others to do so.
+# ========================================================================================
+
+<parthenon/job>
+problem_id = stepdiff
+
+<parthenon/mesh>
+refinement = none
+
+nx1 = 100
+x1min = -0.5
+x1max = 0.5
+ix1_bc = outflow
+ox1_bc = outflow
+
+nx2 = 1
+x2min = -0.5
+x2max = 0.5
+ix2_bc = periodic
+ox2_bc = periodic
+
+nx3 = 1
+x3min = -0.5
+x3max = 0.5
+ix3_bc = periodic
+ox3_bc = periodic
+
+<parthenon/swarm>
+ix1_bc = jaybenne_reflecting
+ox1_bc = jaybenne_reflecting
+ix2_bc = periodic
+ox2_bc = periodic
+ix3_bc = periodic
+ox3_bc = periodic
+
+<parthenon/meshblock>
+nx1 = 100
+nx2 = 1
+nx3 = 1
+
+<parthenon/time>
+tlim = 3.335641e-10
+integrator = rk1
+
+<jaybenne>
+use_ddmc = true
+num_particles = 100000
+dt = 3.335641e-11
+n_nubins = 16
+numin = 1.e12
+numax = 1.e17
+do_emission = false
+do_feedback = false
+transport_model = zone_size
+tracking_algo = history # event
+source_strategy = uniform
+seed = 349857
+
+<mcblock>
+frequency_type = multigroup
+specific_heat = 1.0e8 # erg/K/g
+initial_density = 1.0 # g cm^-3
+initial_temperature = 1.0e5 # K
+initial_radiation = thermal
+
+<mcblock/absorption>
+opacity_model = none
+
+<mcblock/scattering>
+opacity_model = powerlaw
+kappa0 = 0.75e3
+rho_exp = 0.0
+temp_exp = 0.0
+nu_exp = -0.25
+nu_ref = 8.0584e15
+nu_off = 4.5316e14
+rho_ref = 1.0
+temp_ref = 1.0
+
+<parthenon/output0>
+file_type = hdf5
+dt = 3.335641e-10
+variables = field.material.density,         &
+            field.material.sie,             &
+            field.material.internal_energy, &
+            field.jaybenne.energy_tally
diff --git a/inputs/stepdiff_smr.in b/inputs/stepdiff_smr.in
@@ -76,14 +76,18 @@ seed = 349857
 
 <mcblock>
 frequency_type = gray
-opacity_model = none
-scattering_model = constant
-scattering_constant_value = 1.0e3
 specific_heat = 1.0e8 # erg/K/g
 initial_density = 1.0 # g cm^-3
 initial_temperature = 1.0e5 # K
 initial_radiation = thermal
 
+<mcblock/absorption>
+opacity_model = none
+
+<mcblock/scattering>
+opacity_model = constant
+constant_value = 1.0e3
+
 <parthenon/output0>
 file_type = hdf5
 dt = 3.335641e-10

diff --git a/inputs/stepdiff_smr_ddmc.in b/inputs/stepdiff_smr_ddmc.in
@@ -78,14 +78,18 @@ seed = 349857
 
 <mcblock>
 frequency_type = gray
-opacity_model = none
-scattering_model = constant
-scattering_constant_value = 1.0e3
 specific_heat = 1.0e8 # erg/K/g
 initial_density = 1.0 # g cm^-3
 initial_temperature = 1.0e5 # K
 initial_radiation = thermal
 
+<mcblock/absorption>
+opacity_model = none
+
+<mcblock/scattering>
+opacity_model = constant
+constant_value = 1.0e3
+
 <parthenon/output0>
 file_type = hdf5
 dt = 3.335641e-10

diff --git a/inputs/stepdiff_smr_hybrid.in b/inputs/stepdiff_smr_hybrid.in
@@ -78,14 +78,18 @@ diagnostic_level = 1
 
 <mcblock>
 frequency_type = gray
-opacity_model = none
-scattering_model = constant
-scattering_constant_value = 1.0e3
 specific_heat = 1.0e8 # erg/K/g
 initial_density = 1.0 # g cm^-3
 initial_temperature = 1.0e5 # K
 initial_radiation = thermal
 
+<mcblock/absorption>
+opacity_model = none
+
+<mcblock/scattering>
+opacity_model = constant
+constant_value = 1.0e3
+
 <parthenon/output0>
 file_type = hdf5
 dt = 3.335641e-10
+10 −3		README.md
+18 −0		plan_histories/69.md
+25 −0		plan_histories/70.md
+9 −6		singularity-opac/neutrinos/thermal_distributions_neutrinos.hpp
+94 −18		singularity-opac/photons/powerlaw_opacity_photons.hpp
+125 −0		singularity-opac/photons/powerlaw_s_opacity_photons.hpp
+5 −2		singularity-opac/photons/s_opac_photons.hpp
+11 −5		singularity-opac/photons/thermal_distributions_photons.hpp
+3 −1		test/CMakeLists.txt
+32 −29		test/test_mean_opacities.cpp
+239 −1		test/test_powerlaw_opacities.cpp
+2 −2		test/test_spiner_opac_neutrinos.cpp
+85 −1		test/test_thermal_dist.cpp