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

• I have performed a self-review of my own code
• 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)

[gridley]

IMO, this should not be its own type of cell, but rather a functionality introduced to the cell base class. I think so far, the role of a "cell" has been to delineate a region of space, and store some extra values.

There is no reason why a DAGCell should not be able to do chord length sampling, but that functionality has been artificially restricted here. The sub-cell material model might be better off as an empty pointer added to the cell class. There are multiple ways to have a sub-cell heterogeneity model. Exponential distribution chord length sampling is one of them, and I believe there are other stochastic media models out there which would use a similar interface.

I'd love to hear someone else's opinion on how this should be architected. Maybe we should clarify what the abstract role of a Cell is a little better.

[gridley]

These are all unnecessary with public inheritance.

[cn-skywalker]

Thank you very much for your guidance. According to your suggestion, I will rethink how to construct this method.

My current plan is to first implement the CLS method in CSGCell and then expand it to DAGCell after verification.

My preliminary idea is as follows:

1. Add a new abstract class called subcellModel, and it will be added to the cell class as an empty pointer.
2. Create a new class called CLSModel, which inherits from the subcellModel class. All calculations related to CLS chord length sampling will be implemented in this class.
3. Modify the transport logic to implement the CLS method.

If you have any other suggestions, please inform me in a timely manner. I will do my utmost to complete them.

[cn-skywalker]

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

[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?