Fixing double delta_impurity scaling - #56
Conversation
WalkthroughUpdates irradiation scenarios, reaction-rate and activity calculations, near-field dose JSON output, diagnostic logging, regression tests, and documentation for the ChangesCalculation pipeline
Estimated code review effort: 3 (Moderate) | ~25 minutes Suggested reviewers: Poem
🚥 Pre-merge checks | ✅ 4 | ❌ 1❌ Failed checks (1 warning)
✅ Passed checks (4 passed)
✨ Finishing Touches📝 Generate docstrings
🧪 Generate unit tests (beta)
Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out. Comment |
…d dose comparison with RadPro in test
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
There was a problem hiding this comment.
Actionable comments posted: 3
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Inline comments:
In `@docs/source/examples.rst`:
- Line 93: Revise the documentation statements around the “All results” and FM
descriptions to distinguish the default 1 g, per-gram output from values scaled
to the mass supplied with --m. State that mass-scaled FM represents total
photons per second for the requested component mass, and update the
corresponding occurrences near the other referenced examples consistently.
In `@src/f4epurity/dose.py`:
- Around line 263-299: Update the line-source detection and perpendicular
sampling in the dose calculation block: determine whether the source is a line
from the endpoint coordinates (x1, y1, z1) and (x2, y2, z2), not from len(dose),
and construct each test point by applying a normalized offset perpendicular to
the line direction. Preserve the existing point-source formula only when no
valid distinct endpoint pair is provided, and continue storing line results with
the *_cm_perpendicular keys.
In `@tests/test_cobalt_activity.py`:
- Line 95: Remove the redundant f-string prefixes from the three static print
messages in the affected test code, including the message beginning “✓ F4Epurity
calculation complete:” and the corresponding strings at the other referenced
locations, while preserving their output unchanged.
🪄 Autofix (Beta)
Fix all unresolved CodeRabbit comments on this PR:
- Push a commit to this branch (recommended)
- Create a new PR with the fixes
ℹ️ Review info
⚙️ Run configuration
Configuration used: Organization UI
Review profile: CHILL
Plan: Pro
Run ID: 5b206501-2182-44ed-ba58-78eb5ba59484
📒 Files selected for processing (8)
docs/source/examples.rstsrc/f4epurity/decay_chain_calc.pysrc/f4epurity/dose.pysrc/f4epurity/main.pysrc/f4epurity/parsing.pysrc/f4epurity/reaction_rate.pytests/test_cobalt_activity.pytests/test_decay_chain_calc.py
| If multiple sources are provided, by default, it is assumed that all the components where the impurities are located have the same mass. In this way, the dose deviation maps can be easily summed together, and the results are provided as mSv/h/g of the component. If this assumption is not valid, masses (in g) for each component need to be provided through the ``--m`` option. This data can also be provided through the ``--sources_csv`` option by adding a column ``m`` to the `.csv` file. | ||
|
|
||
| **All results are output per gram of material.** No newline at end of file | ||
| **All results are output per gram of material.** |
There was a problem hiding this comment.
🗄️ Data Integrity & Integration | 🟠 Major | ⚡ Quick win
Clarify FM units when --m is supplied.
These statements imply that the generated FM value is always per gram, but main.py passes the requested mass to PointSource, and PointSource._compute_inventory() scales activities before GlobalPointSource.to_sdef() sums them. With --m, FM is therefore the total photons/second for that component mass; only the default 1 g case is numerically per gram.
Please revise the “all results” and FM descriptions to distinguish the default 1 g output from mass-scaled output.
Also applies to: 104-112, 140-140
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In `@docs/source/examples.rst` at line 93, Revise the documentation statements
around the “All results” and FM descriptions to distinguish the default 1 g,
per-gram output from values scaled to the mass supplied with --m. State that
mass-scaled FM represents total photons per second for the requested component
mass, and update the corresponding occurrences near the other referenced
examples consistently.
| if hasattr(dose, "__iter__") and len(dose) > 1: | ||
| # Line source case - calculate at perpendicular distances | ||
| logging.info("Line Source Mode:") | ||
| line_center_x = (x1 + x2) / 2 | ||
| line_center_y = (y1 + y2) / 2 | ||
| line_center_z = (z1 + z2) / 2 | ||
|
|
||
| for dist in specific_distances: | ||
| # Calculate dose at perpendicular distance from line center | ||
| # Move in x-direction for simplicity (could be made more general) | ||
| test_point_x = line_center_x + dist | ||
| dose_at_distance = dose_from_line_source( | ||
| dose, x1, y1, z1, x2, y2, z2, | ||
| test_point_x, line_center_y, line_center_z | ||
| ) | ||
| dose_at_distances[f"{dist}_cm_perpendicular"] = float(dose_at_distance) | ||
| logging.info(f"Perpendicular distance: {dist:6.1f} cm -> Dose: {dose_at_distance:.3e} μSv/h/g") | ||
| else: | ||
| # Point source case | ||
| logging.info("Point Source Mode:") | ||
| for dist in specific_distances: | ||
| # Point source: use 1/r^2 formula | ||
| dose_at_distance = dose[0] / (4 * pi * dist**2) | ||
| dose_at_distances[f"{dist}_cm"] = float(dose_at_distance) | ||
| logging.info(f"Distance: {dist:6.1f} cm -> Dose: {dose_at_distance:.3e} μSv/h/g") | ||
|
|
||
|
|
||
|
|
||
| # Save dose at specific distances to a JSON file | ||
| if x2 is not None and y2 is not None and z2 is not None: | ||
| dose_distances_file = f"{run_dir}/dose_at_distances_{x1}_{y1}_{z1}_to_{x2}_{y2}_{z2}.json" | ||
| else: | ||
| dose_distances_file = f"{run_dir}/dose_at_distances_{x1}_{y1}_{z1}.json" | ||
| with open(dose_distances_file, "w") as f: | ||
| json.dump(dose_at_distances, f, indent=4) | ||
| logging.info(f"Saved dose at specific distances to: {dose_distances_file}") | ||
|
|
There was a problem hiding this comment.
🎯 Functional Correctness | 🟠 Major | ⚡ Quick win
Derive a true perpendicular offset and identify line sources from their endpoints.
line_center_x + dist is only perpendicular for a line with no x component. For an x-aligned line it remains on the source, so the reported *_cm_perpendicular dose is wrong. Also, a one-cell line source enters the point-source branch because detection depends on len(dose).
Proposed fix
- if hasattr(dose, "__iter__") and len(dose) > 1:
+ is_line_source = all(value is not None for value in (x2, y2, z2))
+ if is_line_source:
logging.info("Line Source Mode:")
line_center_x = (x1 + x2) / 2
line_center_y = (y1 + y2) / 2
line_center_z = (z1 + z2) / 2
+ direction = np.array([x2 - x1, y2 - y1, z2 - z1], dtype=float)
+ reference_axis = np.array([1.0, 0.0, 0.0])
+ if np.allclose(np.cross(direction, reference_axis), 0):
+ reference_axis = np.array([0.0, 1.0, 0.0])
+ perpendicular = np.cross(direction, reference_axis)
+ perpendicular /= np.linalg.norm(perpendicular)
for dist in specific_distances:
- test_point_x = line_center_x + dist
+ test_point_x = line_center_x + dist * perpendicular[0]
+ test_point_y = line_center_y + dist * perpendicular[1]
+ test_point_z = line_center_z + dist * perpendicular[2]
dose_at_distance = dose_from_line_source(
dose, x1, y1, z1, x2, y2, z2,
- test_point_x, line_center_y, line_center_z
+ test_point_x, test_point_y, test_point_z
)📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
| if hasattr(dose, "__iter__") and len(dose) > 1: | |
| # Line source case - calculate at perpendicular distances | |
| logging.info("Line Source Mode:") | |
| line_center_x = (x1 + x2) / 2 | |
| line_center_y = (y1 + y2) / 2 | |
| line_center_z = (z1 + z2) / 2 | |
| for dist in specific_distances: | |
| # Calculate dose at perpendicular distance from line center | |
| # Move in x-direction for simplicity (could be made more general) | |
| test_point_x = line_center_x + dist | |
| dose_at_distance = dose_from_line_source( | |
| dose, x1, y1, z1, x2, y2, z2, | |
| test_point_x, line_center_y, line_center_z | |
| ) | |
| dose_at_distances[f"{dist}_cm_perpendicular"] = float(dose_at_distance) | |
| logging.info(f"Perpendicular distance: {dist:6.1f} cm -> Dose: {dose_at_distance:.3e} μSv/h/g") | |
| else: | |
| # Point source case | |
| logging.info("Point Source Mode:") | |
| for dist in specific_distances: | |
| # Point source: use 1/r^2 formula | |
| dose_at_distance = dose[0] / (4 * pi * dist**2) | |
| dose_at_distances[f"{dist}_cm"] = float(dose_at_distance) | |
| logging.info(f"Distance: {dist:6.1f} cm -> Dose: {dose_at_distance:.3e} μSv/h/g") | |
| # Save dose at specific distances to a JSON file | |
| if x2 is not None and y2 is not None and z2 is not None: | |
| dose_distances_file = f"{run_dir}/dose_at_distances_{x1}_{y1}_{z1}_to_{x2}_{y2}_{z2}.json" | |
| else: | |
| dose_distances_file = f"{run_dir}/dose_at_distances_{x1}_{y1}_{z1}.json" | |
| with open(dose_distances_file, "w") as f: | |
| json.dump(dose_at_distances, f, indent=4) | |
| logging.info(f"Saved dose at specific distances to: {dose_distances_file}") | |
| is_line_source = all(value is not None for value in (x2, y2, z2)) | |
| if is_line_source: | |
| logging.info("Line Source Mode:") | |
| line_center_x = (x1 + x2) / 2 | |
| line_center_y = (y1 + y2) / 2 | |
| line_center_z = (z1 + z2) / 2 | |
| direction = np.array([x2 - x1, y2 - y1, z2 - z1], dtype=float) | |
| reference_axis = np.array([1.0, 0.0, 0.0]) | |
| if np.allclose(np.cross(direction, reference_axis), 0): | |
| reference_axis = np.array([0.0, 1.0, 0.0]) | |
| perpendicular = np.cross(direction, reference_axis) | |
| perpendicular /= np.linalg.norm(perpendicular) | |
| for dist in specific_distances: | |
| # Calculate dose at perpendicular distance from line center | |
| # Move in x-direction for simplicity (could be made more general) | |
| test_point_x = line_center_x + dist * perpendicular[0] | |
| test_point_y = line_center_y + dist * perpendicular[1] | |
| test_point_z = line_center_z + dist * perpendicular[2] | |
| dose_at_distance = dose_from_line_source( | |
| dose, x1, y1, z1, x2, y2, z2, | |
| test_point_x, test_point_y, test_point_z | |
| ) | |
| dose_at_distances[f"{dist}_cm_perpendicular"] = float(dose_at_distance) | |
| logging.info(f"Perpendicular distance: {dist:6.1f} cm -> Dose: {dose_at_distance:.3e} μSv/h/g") | |
| else: | |
| # Point source case | |
| logging.info("Point Source Mode:") | |
| for dist in specific_distances: | |
| # Point source: use 1/r^2 formula | |
| dose_at_distance = dose[0] / (4 * pi * dist**2) | |
| dose_at_distances[f"{dist}_cm"] = float(dose_at_distance) | |
| logging.info(f"Distance: {dist:6.1f} cm -> Dose: {dose_at_distance:.3e} μSv/h/g") | |
| # Save dose at specific distances to a JSON file | |
| if x2 is not None and y2 is not None and z2 is not None: | |
| dose_distances_file = f"{run_dir}/dose_at_distances_{x1}_{y1}_{z1}_to_{x2}_{y2}_{z2}.json" | |
| else: | |
| dose_distances_file = f"{run_dir}/dose_at_distances_{x1}_{y1}_{z1}.json" | |
| with open(dose_distances_file, "w") as f: | |
| json.dump(dose_at_distances, f, indent=4) | |
| logging.info(f"Saved dose at specific distances to: {dose_distances_file}") |
🧰 Tools
🪛 ast-grep (0.44.1)
[warning] 295-295: File path is request-/variable-derived; validate and normalize to prevent path traversal.
Context: open(dose_distances_file, "w")
Note: [CWE-22] Improper Limitation of a Pathname to a Restricted Directory ('Path Traversal').
(open-filename-from-request)
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In `@src/f4epurity/dose.py` around lines 263 - 299, Update the line-source
detection and perpendicular sampling in the dose calculation block: determine
whether the source is a line from the endpoint coordinates (x1, y1, z1) and (x2,
y2, z2), not from len(dose), and construct each test point by applying a
normalized offset perpendicular to the line direction. Preserve the existing
point-source formula only when no valid distinct endpoint pair is provided, and
continue storing line results with the *_cm_perpendicular keys.
| tolerance = 0.01 # 1% | ||
|
|
||
| # Print result for user visibility | ||
| print(f"\n✓ F4Epurity calculation complete:") |
There was a problem hiding this comment.
📐 Maintainability & Code Quality | 🟡 Minor | ⚡ Quick win
Remove redundant f-string prefixes.
These strings contain no interpolation; Ruff reports F541, which can fail linting.
Proposed fix
- print(f"\n✓ F4Epurity calculation complete:")
+ print("\n✓ F4Epurity calculation complete:")
...
- pytest.fail(f"Could not find dose conversion factor for Co-60 in any naming format")
+ pytest.fail("Could not find dose conversion factor for Co-60 in any naming format")
...
- print(f"\n✓ Dose calculation results:")
+ print("\n✓ Dose calculation results:")Also applies to: 123-123, 147-147
🧰 Tools
🪛 Ruff (0.15.21)
[error] 95-95: f-string without any placeholders
Remove extraneous f prefix
(F541)
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In `@tests/test_cobalt_activity.py` at line 95, Remove the redundant f-string
prefixes from the three static print messages in the affected test code,
including the message beginning “✓ F4Epurity calculation complete:” and the
corresponding strings at the other referenced locations, while preserving their
output unchanged.
Source: Linters/SAST tools
|
This pull request wants to close this issue #55 |
Radiation-Transport
assigned
Radiation-Transport and
alexvalentine94 and unassigned
Radiation-Transport
|
Thanks @digiacomo4 & @Radiation-Transport. I will look into this no later than 29.07. @teade for information. |
There was a problem hiding this comment.
Thanks Matteo for this work. Having looked at this independently I would agree that there was a bug in double counting the delta_impurity. Reviewing the report, I can see that while we did relatively rigorous V&V they were focused on specific aspects of the code rather than complete end to end tests which you have integrated.
@teade also performed an independent check and agrees.
I have left some minor comments for consideration.
|
|
||
| - The source position (x, y, z coordinates) | ||
| - The photon energy spectrum from the activated impurities | ||
| - The photon emission rate (FM value) expressed in photons per second per gram of material |
There was a problem hiding this comment.
Be more explicit here, 'first entry on the FM, tally multiplier card, in MCNP'
| # REF : ITER_D_8WK64Y | ||
| # old SA2 and generic test and 1 year of irradiation scenarios included for testing purposes. | ||
| # REF for DT1 and DT2: CXM7AR v1.3 Safety irradiation scenario | ||
| IRRAD_SCENARIOS = { |
There was a problem hiding this comment.
Have not checked the DT1 and DT2 scenarios. Will assume correct but please double check as these are hard coded.
|
|
||
| # Calculate the reaction rate | ||
| reaction_rate = (delta_impurity / 100) * sigma * total_flux * 1e-24 | ||
| reaction_rate = sigma * total_flux * 1e-24 |
There was a problem hiding this comment.
Agree that we were double counting and that this should be removed from here.
|
|
||
| # Calculate the reaction rate | ||
| reaction_rate = (delta_impurity / 100) * sigma_eff * total_flux * 1e-24 | ||
| reaction_rate = sigma_eff * total_flux * 1e-24 |
| sigma_value = float(sigma_eff.flat[0]) | ||
| else: | ||
| sigma_value = float(sigma_eff) | ||
| logging.info(f" Reaction {parent} -> {product}: sigma_eff = {sigma_value:.6e} barns") |
There was a problem hiding this comment.
Most of the effective cross sections we looked at in the report are 0.1 -> 100. Is the 6 decimal place printing excessive?
| logging.info(f" Flux spectrum: {flux_spectrum}") | ||
|
|
||
| # Calculate the reaction rate based on the flux and effective cross section | ||
| reaction_rate = calculate_reaction_rate( |
There was a problem hiding this comment.
we can remove delta_impurity from the function.
| activity = ( | ||
| final_nuclides[nuclide] * decay_data_dic[nuclide]["lambda"][0] | ||
| ) | ||
| n_atoms = final_nuclides[nuclide] |
There was a problem hiding this comment.
OK. Just readability changes and addition of logging.
| specific_distances = [1, 10, 50, 100, 200] # distances in cm | ||
| dose_at_distances = {} | ||
|
|
||
| logging.info("\n--- Dose at Specific Distances ---") |
There was a problem hiding this comment.
Looks like this is output by default, do we want to make optional? Happy either way.
Perhaps we should also add a couple of lines in the documentation on this feature and why it has been added.
| reaction_rate_co60 = calculate_reaction_rate(delta_impurity, sigma_eff_co60, flux_spectrum) | ||
| reaction_rate_co60m = calculate_reaction_rate(delta_impurity, sigma_eff_co60m, flux_spectrum) | ||
|
|
||
| # Calculate number of atoms for Co-59 |
There was a problem hiding this comment.
Do this with F4Epurity rather than manually?
There was a problem hiding this comment.
Thanks for this - it is very good to have these integration tests. My only comment is that all of the print statements are unnecessary. With the unit/integration testing we are really only concerned if the test passes/fails. These are just additional noise when you run pytest.
| sigma_eff_co60m = 3.873241 # barns | ||
|
|
||
| # Calculate reaction rates using F4Epurity function | ||
| reaction_rate_co60 = calculate_reaction_rate(delta_impurity, sigma_eff_co60, flux_spectrum) |
There was a problem hiding this comment.
Can remove the delta_impurity variable
|
|
||
| # Calculate reaction rates using F4Epurity function | ||
| reaction_rate_co60 = calculate_reaction_rate(delta_impurity, sigma_eff_co60, flux_spectrum) | ||
| reaction_rate_co60m = calculate_reaction_rate(delta_impurity, sigma_eff_co60m, flux_spectrum) |
|
Once approved, we can bump to v1.1.0 |
3 participants
Summary by CodeRabbit
New Features
--write_sdef, including source position, photon spectrum, emission rates, mass defaults, and FM scaling.Bug Fixes
Tests