Ability to compute material volume fractions over mesh elements

[paulromano]

Description

This PR introduces a capability to compute the volume of individual materials within a mesh element. The intended purpose of this feature is to support mesh-based R2S workflows. Namely, this capability can be used to compute homogenized material compositions over mesh elements, which are then activated to create a decay photon source over a mesh. A few notes:

• I was able to reuse some of the logic from volume calculations by separating out a reduce_indices_hits function that is now used both for regular volume calculations as well as the new feature here.
• To support this feature, two new C API functions are added: openmc_mesh_get_n_elements and openmc_mesh_volume_fractions.
• Users can utilize this functionality through openmc.lib as the Mesh.volume_fractions(...) method.
• There was an indexing bug in RectilinearMesh::volume that is fixed here.
• Void volumes are calculated as well and are represented with None on the Python side

One of the challenging design aspects of this was how to handle memory allocation, since we don't know a priori how many materials are going to be hit in a given mesh element. The way I chose to handle this is as follows. First, the loop over mesh elements is done on the Python side. This avoids the need to handle a doubly nested allocation (akin to <vector<vector<MaterialVolume>>) and instead means we only need a single allocation per call to openmc_mesh_volume_fractions. This allocation is done on the Python side with ctypes with a guess on how many materials at most will be present (starts at 16). If you really run into a situation where a single element has more than 16 materials, it will raise an exception on the Python side, which then triggers a reallocation with double the size.

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)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[pshriwise]

Pleasantly surprised at the small amount of code this takes. Thanks @paulromano!! Not much to clean up here.

One design comment/thought that came to mind was that it might be useful to wrap the current version of Mesh::volume_fractions around a function that returns a vector of MaterialVolumes so we can suport return of dynamically sized results to support cases where this capability may be accessed via the C++ API.

[paulromano]

@pshriwise All your feedback has been addressed. Note that I ended up changing the name of the method to Mesh::material_volumes since the method returns volumes, not volume fractions (returning volumes gives more information, and volume fractions can be calculated trivially after the fact). I also added a (currently unused) overload for Mesh::material_volumes that returns a vector as you had requested.

[shimwell]

Perhaps I've misunderstood something here but I tried to make a minimal example and was unsure about the result.

I thought the material mass would be 1/8 of the volume of the sphere as my mesh with 8 voxels covers the sphere, but I'm getting 1.0 for each material volume.

import math
import openmc
import openmc.lib
mat_1 = openmc.Material(material_id=1) 
mat_1.add_element('Ag', 1)
mat_1.set_density('g/cm3', 1)
mat_2 = openmc.Material(material_id=2) 
mat_2.add_element('Au', 1)
mat_2.set_density('g/cm3', 1)
sphere_radius = 1
materials = openmc.Materials([mat_1, mat_2])
sph1 = openmc.Sphere(r=sphere_radius, boundary_type='vacuum')
zplane1 = openmc.ZPlane(z0=0)
reg1 = -sph1 & +zplane1
reg2 = -sph1 & -zplane1
cell1 = openmc.Cell(region=reg1)
cell2 = openmc.Cell(region=reg2)
cell1.fill = mat_1
cell2.fill = mat_2
geometry = openmc.Geometry([cell1, cell2])
settings = openmc.Settings()
settings.batches = 1
settings.particles = 1
settings.run_mode = 'fixed source'
model = openmc.Model(geometry=geometry, materials=materials, settings=settings)
model.export_to_model_xml()
openmc.lib.init()
mesh = openmc.lib.RegularMesh()
mesh.dimension = (2, 2, 2)  # 8 elements
mesh.set_parameters(
    lower_left=(-1, -1, -1), # fully covers the sphere of radius 1
    upper_right=(1, 1, 1),
)
volume_of_sphere = (4/3) * math.pi * math.pow(sphere_radius, 3)
vols = mesh.material_volumes(n_samples = 100000)
for vol in vols:
    print(vol)
    # assert vol[0][1] == volume_of_sphere/mesh.n_elements

[paulromano]

@shimwell Your example there helped to uncover a limitation, which was that we had implicitly assumed that the mesh was fully contained within the defined geometry. That limitation has now been removed. While looking at the surrounding code, I also realized that we were not handle random number seeds correctly in parallel, so that has also been fixed.

[shimwell]

Thanks for the fix Paul. Running the above script on the PR branch is now producing expected results.

I see there are tests in ```tests/unit_tests/test_lib.py`` that check the material volumes for a few situations.

I can see a tests where the material volumes are 100% of the mesh cell.
I can see a tests where there is a single material in the mesh voxel.

Wondering if it is worth adding a test where there are 2 materials and they don't fill the voxel 100%

[paulromano]

@shimwell This snippet covers the case you just mentioned (multiple materials that don't fill a voxel 100%):

https://github.com/paulromano/openmc/blob/0edd9ba29b82b088ad84f0df461e8e2e7b0b0a16/tests/unit_tests/test_lib.py#L594-L601

[shimwell]

This looks good to me. I'm now wondering how you make use of this in the r2S workflow. Perhaps best if @pshriwise approves this PR

[paulromano]

I'll have some mesh-based R2S examples to share in the near-term that use this capability.

@pshriwise Any further requests on this PR?

[pshriwise]

No further comments from me! I'll give this until tomorrow morning as it's a substantial feature addition. Looking forward to trying it out @paulromano! Thanks!