New Feature Development: The Implementation of chord length sampling (CLS) method for stochastic media( #3286 )

[cn-skywalker]

Description

In #3286, the Chord Length Sampling (CLS) method was described for enhancing modeling capabilities in OpenMC. This PR begins the implementation of this feature.

The main changes in this PR include:

1. Modifications to the Cell class: Added particle_radius_ and fill_ratio_ attributes to store the radius and fill ratio of dispersed particles.
2. New CLSCell class: Inherits from the CSGCell class. Added the set_translation function to dynamically modify translation_.
3. Modified read_cells function: Updated to initialize the CLSCell class.

Test Case

This test case describes a scenario where Universe 1 (composed of cell1 and cell2) is randomly dispersed in cell3. A new dispersion tag is added to store the particle_radius and fill_ratio variables.

<?xml version='1.0' encoding='UTF-8'?>
<geometry>
  <cell id="1" material="3" name="sphere cell" region="-1" universe="1"/>
  <cell id="2" material="1" name="base cell" region="1" universe="1"/>
  <cell fill="1" id="3" name="fuel cell" region="2 -3 4 -5" translation="0.1 0.1 0.0" universe="2">
      <dispersion particle_radius="0.05" fill_ratio="0.3"/>
  </cell>
  <cell id="4" material="2" name="water cell" region="~(2 -3 4 -5) (6 -7 8 -9)" universe="2"/>
  <surface coeffs="1.0 1.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 -0.01" id="1" name="sphere_surface" type="quadric"/>
  <surface coeffs="-0.5" id="2" name="minimum x" type="x-plane"/>
  <surface coeffs="0.5" id="3" name="maximum x" type="x-plane"/>
  <surface coeffs="-0.5" id="4" name="minimum y" type="y-plane"/>
  <surface coeffs="0.5" id="5" name="maximum y" type="y-plane"/>
  <surface boundary="reflective" coeffs="-1.0" id="6" name="minimum x" type="x-plane"/>
  <surface boundary="reflective" coeffs="1.0" id="7" name="maximum x" type="x-plane"/>
  <surface boundary="reflective" coeffs="-1.0" id="8" name="minimum y" type="y-plane"/>
  <surface boundary="reflective" coeffs="1.0" id="9" name="maximum y" type="y-plane"/>
</geometry>

Test Results

Testing confirmed that the newly added CLSCell class can correctly read the dispersion tag.

Future Development Plan

To minimize modifications to the existing transport logic, I plan to implement on-the-fly modeling of dispersed particles by dynamically adjusting the translation_ property of cell3 during runtime.

Acknowledgments

As this is my first attempt at developing such a significant feature, any discussions and guidance would be greatly appreciated!

Checklist

• [x] I have performed a self-review of my own code
• [x] I have run clang-format (version 15) on any C++ source files (if applicable)
• [ ] I have followed the style guidelines for Python source files (if applicable)

[pshriwise]

Thanks for this contribution @cn-skywalker! It's exciting that we can make this method available. It's come across my desk a few times and I like it.

I'll second @gridley's comments on generalizing the method to support both CSG and DAGMC/CAD cell types. The additional pointer specifying that the cell is filled with a CLS distribution could be a good approach. I'm thinking of this as an alternative cell fill type in the hope that we can implement it in this way without affecting higher level geometry methods (i.e. find_cell). Looking forward to talking through it!

By the way, do I need to close this PR and start another one, or should I continue to make modifications in this PR?

Please feel free to make updates in the same PR, yes! What I'll do for now is make this a draft PR that you can mark as ready for review when the time comes.

[cn-skywalker]

Currently, I have developed the Stochastic_Media class modeled after the Material class, and plan to implement it within the material input card section.I performed a local reset to the develop branch to revert previous inappropriate modifications, resulting in a non-linear commit history. The latest three commits currently reside on the develop branch as functional updates.

[cn-skywalker]

Description

I added test input cards for CLS_media which demonstrate that the relevant parameters can be correctly read. The next step will be to implement stochastic media as a class similar to materials and universes, which can be filled into cells.

Test Case Example

  <CLS_media id="400" particle_material="100" matrix_material="200" radius="0.1" pack_fraction="0.3"/>
  <CLS_media id="500" particle_material="100" matrix_material="200" radius="0.1" pack_fraction="0.3"/>

Next Plan

Enable stochastic media to be assigned to cells like materials and universes

[cn-skywalker]

Development of Stochastic Media Functionality

(1) Integration as Cell Filling Type

• Treat stochastic media as a cell filling type alongside existing types (material, universe, and lattice).

(2) Modifications to find_cell_inner Function

• When particles transport into cells filled with stochastic media:
  • Randomly select material composition based on filling ratios defined in the stochastic media
  • Update particle's status_ parameter to flag stochastic media interaction
  • This flag will be used for stochastic media sampling in event_advance()

(3) Modifications to event_advance Function

• When particle status indicates stochastic media interaction:
  • Add stochastic media distance sampling
  • If this distance is the minimum: Perform a virtual surface crossing in transport_history_based_single_particle Function

(4) Modifications to transport_history_based_single_particle Function

• When particles are flagged for virtual surface crossing:
  • Update particle material based on stochastic media information
  • Current implementation: Basic functionality supporting only:
    • One type of particulate material
    • One type of matrix material

[cn-skywalker]

3. Verification Cases

I have added the following test cases for validation:

(1) Comparison with Homogenization Method (VHM)

• Created a single fuel pin model:
  • VHM model: Fuel cell fill with UO₂ material
  • CLS model:
    • Particulate material: UO₂ (radius=0.02cm)
    • Matrix material: UO₂
    • Filling ratio: 10%
• Results: VHM: kinf = 1.62818 CLS: kinf = 1.62825

This result looks great, but the following validation confuses me

(2) Comparison with Explicit Particulate Model (RSA)

• Created benchmark model using RSA method with identical parameters:
• Matrix material: Zirconium
• Particulate material: UO₂
• Same geometry and filling ratio as CLS model
• Results: RSA (reference): kinf = 1.61606 CLS: kinf = 1.64832

I've put the corresponding input cards in this zip file model.zip

❓ Unexpected Behavior

The significant discrepancy between CLS and RSA results raises several concerns:

1. The CLS result shows higher kinf than VHM when using Zirconium matrix, contrary to expectations. This phenomenon does not make sense. Even if the string length is not reasonable, it does not explain why the kinf became higher after the matrix changed from UO2 to zirconium. I find this inexplicable.
2. Detailed particle tracking during debugging showed:

• No anomalies in transport behavior
• Correct material sampling in stochastic media
• Proper virtual surface crossing implementation

🆘 Community Assistance Request

@pshriwise @gridley This question has puzzled me for a long time. Any community discussions or assistance would be greatly appreciated! 🙏

[gridley]

Very interesting, thank you for sharing your detailed findings! Maybe running a stochastic volume calculation to check that the material volumes are truly equivalent would be one step towards a solid match against the explicit geometry model. It seems that the difference in k is larger than you'd typically expect from a fine-scale heterogeneity effect, so I would guess this relates to material volume conservation.

I'd love to dive deeper into this method with you in the future, but won't be able to help much on this for now, sorry.

[cn-skywalker]

Maybe running a stochastic volume calculation to check that the material volumes are truly equivalent would be one step towards a solid match against the explicit geometry model. It seems that the difference in k is larger than you'd typically expect from a fine-scale heterogeneity effect, so I would guess this relates to material volume conservation.

Additional Debugging Details

Volume Analysis for UO₂ Dispersed in Zirconium

• Material Definitions:
  • Material 11: Dispersed particulate material (UO₂)
  • Material 12: Matrix material (Zirconium)

Results Comparison

Method	Material 11 Volume (cm³)	Material 12 Volume (cm³)	k Result
RSA	0.04215680499575621	0.3794961144660949	1.61606
CLS	0.042165008431963885	0.37948791102988716	1.64832

Visualization

• Figure 1: RSA Benchmark Results
  

• Figure 2: CLS Method Results
  

---

@pshriwise @gridley The results showed that the difference in material volume among the methods was not significant.

Now only modified the material switching logic in related functions .No changes made to other core modules. Have I overlooked any critical operations during this modification?Any other suggestions to help locate the bug?

[gridley]

Hm, very strange! I will try building a TRISO model and give this a shot when I get the time, although I will warn that it might be a while.

If you want to continue debugging, I would suggest adding some quite granular tallies like a fine energy spectrum filter and running with a very large number of particles. If you are able to discern differences aside from just k, it might help you pin down the source of the discrepancy.

Alternatively, maybe outputting chord lengths from a direct geometry calculation to verify the models are truly equivalent would be another approach. And maybe comparing the collision and track length estimates of k, which should be available in the state point file, would also be useful.

[cn-skywalker]

@gridley Thank you very much for your attention to my work. I seem to have identified the source of the bug. I have statistically analyzed the total chord lengths in both the particles and matrix , but the observed ratio does not match theoretical expectations despite the formulas appearing correct. Further analysis is underway. I sincerely appreciate your assistance!

[cn-skywalker]

The previous significant deviation in the CLS method results has been resolved. Due to my unfamiliarity with C++ operator precedence, the particle transport chord length was miscalculated. After corrections, the CLS test case yields k = 1.61925, showing a +300 pcm deviation from the RSA reference solution (1.61606), which aligns with the expected error margin.