Ray tracing volume calculation and geometry testing

[vitmog]

Description

Cost of the current volume fraction calculation via the rejection estimator principally depends on (and inversely proportional to) the target region volume fraction of the bounding box, therefore it provides troubles with computation of small volume fractions. In this PR, it is implemented an alternate estimator based on combination of the Hit-and-Run algorithm [1] and the ray tracing which uniformly covers entire location-direction model domain. The efficiency of this algorithm less depends on the volume fraction and it requires running geometry routines similarly to the neutral particle transport case. That all allow to pursue several goals:

• Efficient estimation of small volume fractions;
• Fast geometry input diagnostics on preliminary stages;
• Efficient geometry routines reliability testing and their performance measurement.

The switch between both the schemes is provided.

By the way, this PR also:

• Introduces the max_iterations parameter what partially fixes #1702;
• Fixes an error with the vol_calc regression test input which was became explicit after the volume estimator's change;
• Brings the results text output format closer to the OpenMC style.

Usage example

It is considered below a simple test model of spheres within spheres containing 7 cells which $20 \times 20 \times 20$ cm bounding box's volume fractions are in range from 5.23e-1 to 5.23e-10.

The standard hitting volume calculation requires either pass estimator_type = 'hit' or default empty:

cell_vol = openmc.VolumeCalculation(cells_list, 100000000, lower_left, upper_right)

It provides the following output:

Total sample size    = 1.0000e+08 
Running cost         = 1.3692e+02 thread-sec
Cost of hitting      = 1.3692e-06 thread-sec/hit

    Cell     1: 0.0000e+00 +/- 0.0000e+00 cm^3   // 4.1888e-06 | R < 0.01 cm sphere
    Cell     2: 4.8800e-03 +/- 6.2482e-04 cm^3   // 4.1888e-03 | R < 0.1 cm sphere in the center
    Cell     3: 4.2010e+00 +/- 1.8328e-02 cm^3   // 4.1846e+00 | 0.1 < R < 1. cm layer - sm.sph.
    Cell     4: 1.3120e-01 +/- 3.2397e-03 cm^3   // 1.2692e-01 | 1. < R < 1.01 cm layer
    Cell     5: 3.2000e-04 +/- 1.6000e-04 cm^3   // 1.2820e-03 | 1.01 < R < 1.0101 cm layer
    Cell     6: 4.1852e+03 +/- 3.9957e-01 cm^3   // 4.1845e+03 | 1.0101 < R < 10 cm layer
    Cell     7: 3.6800e-03 +/- 5.4259e-04 cm^3   // 4.1888e-03 | R < 0.1 cm sphere on periphery

Ray tracing volume calculation requires definition of estimator_type = 'ray' argument:

cell_vol_rt = openmc.VolumeCalculation(cells_list, 40000000, lower_left, upper_right, 'ray')

It provides the following result which has almost the same running cost as the previous one:

Total sample size    = 4.0000e+07 
Running cost         = 1.3440e+02 thread-sec
Rays traced          = 3.0638e+07 
Average ray length   = 1.0764e+00 segments/ray
Cost of tracing      = 4.0753e-06 thread-sec/segment

    Cell     1: 4.1083e-06 +/- 7.3968e-07 cm^3   // 4.1888e-06 | R < 0.01 cm sphere
    Cell     2: 4.0978e-03 +/- 7.1491e-05 cm^3   // 4.1888e-03 | R < 0.1 cm sphere in the center
    Cell     3: 4.1804e+00 +/- 7.1841e-03 cm^3   // 4.1846e+00 | 0.1 < R < 1. cm layer - sm.sph.
    Cell     4: 1.2719e-01 +/- 2.6854e-04 cm^3   // 1.2692e-01 | 1. < R < 1.01 cm layer
    Cell     5: 1.2819e-03 +/- 3.4497e-06 cm^3   // 1.2820e-03 | 1.01 < R < 1.0101 cm layer
    Cell     6: 4.1848e+03 +/- 4.1176e-01 cm^3   // 4.1845e+03 | 1.0101 < R < 10 cm layer
    Cell     7: 4.1990e-03 +/- 9.2531e-05 cm^3   // 4.1888e-03 | R < 0.1 cm sphere on periphery

In this example, the computational gain is obtained on almost all cells but there is no general rule about higher ray tracing efficiency, it depends on configuration and shape of model.

Further works

If needed, this implementation can be relatively easily extended onto the following functionality:

• Surface areas calculation using normals;
• Automated bounding box determination during calculation or complete refuse of its usage;
• Advanced geometry input diagnostics via ray tracing (if known).

Also, as a quite natural way for volume and areas calculation development could be some unification of the mode RunMode::VOLUME with RunMode::FIXED_SOURCE and RunMode::EIGENVALUE for using Filter and, maybe, Tally classes. For example, it could automatically to provide fixing #2367 and #114.

References

1. R.L. Smith. Efficient Monte Carlo Procedures for Generating Points Uniformly Distributed Over Bounded Regions. 1984, Operations Research, 32(6), 1296–1308. https://doi.org/10.1287/opre.32.6.1296

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)
• [x] I have followed the style guidelines for Python source files (if applicable)
• [x] I have made corresponding changes to the documentation (if applicable)
• [x] I have added tests that prove my fix is effective or that my feature works (if applicable)