Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Ability to compute material volume fractions over mesh elements #2802

Merged
merged 18 commits into from
Jan 19, 2024

Conversation

paulromano
Copy link
Contributor

@paulromano paulromano commented Dec 12, 2023

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)

Copy link
Contributor

@pshriwise pshriwise left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

include/openmc/volume_calc.h Outdated Show resolved Hide resolved
include/openmc/volume_calc.h Outdated Show resolved Hide resolved
src/mesh.cpp Outdated Show resolved Hide resolved
src/mesh.cpp Show resolved Hide resolved
src/mesh.cpp Show resolved Hide resolved
src/mesh.cpp Show resolved Hide resolved
@paulromano
Copy link
Contributor Author

@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
Copy link
Member

shimwell commented Jan 4, 2024

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
Copy link
Contributor Author

@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
Copy link
Member

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
Copy link
Contributor Author

@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
Copy link
Member

shimwell commented Jan 17, 2024

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
Copy link
Contributor Author

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?

Copy link
Contributor

@pshriwise pshriwise left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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!

@pshriwise pshriwise added the Merging Soon PR will be merged in < 24 hrs if no further comments are made. label Jan 18, 2024
@pshriwise pshriwise merged commit 187160d into openmc-dev:develop Jan 19, 2024
18 checks passed
@paulromano paulromano deleted the mesh-mat-volume-fracs branch January 21, 2024 17:29
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Merging Soon PR will be merged in < 24 hrs if no further comments are made.
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants