Specify temperature from a field (structured mesh only)

[JoffreyDorville]

Description

This PR introduces a new way to specify temperature using a field based on a spatial mesh. The motivation for this work is to simplify data transfer (and workflow) when working with multi-physics framework (e.g., Cardinal). It is limited to structured mesh now as a way to get started, but the objective is to add more complexity (support to unstructured mesh via the XDG library) with follow-up PRs.

This work is still in progress (draft PR).

I am planning to:

• refine the Python API based on other needs I have to support a generic Field class,
• add more consistency checks in the Field class,
• add unit tests for all structured mesh types,
• provide verification results here.

Do not hesitate to start a conversation here if you are interested in this feature and want to share any ideas!

More details on the implementation

This feature relies on a minimal interface with the mesh class (via the methods get_bin() and distance_to_next_boundary()) which is expected to simplify the connection with unstructured mesh classes once XDG functionalities are available in OpenMC.

In this PR, the temperature field is defined as the association of a spatial mesh and a list of temperature values corresponding to each mesh cell. The temperature in each mesh cell is considered constant in the entire volume of the mesh cell. Now, every time a particle is advancing, OpenMC also checks the distance to the next boundary associated with the mesh of the temperature field. The temperature associated with the particle is updated in three different scenarios (first time the particle is in the transport loop, after crossing a geometric cell boundary, and after crossing a boundary from the mesh of the temperature field) and will take precedence over any temperature value defined for the cell, material, or globally.

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)

[paulromano]

A few initial thoughts:

• I see that this currently only works in history-based mode. You'll need to extend event-based mode as well (see process_advance_particle_events in event.cpp)
• It looks like next_event could just be returned from event_advance() rather than stored as a new attribute on the particle.
• @pshriwise talked about the notion that we may want mesh boundary crossings to be a generic concept (there may be other reasons to stop particles at mesh boundaries). If we were to generalize the event handling, it would involve something like a vector of meshes against which we would check distances to and if any of them are the closest, the "cross mesh" event would be activated. I'm not sure it's worth it to do that as part of this PR though, or even what other meshes we would potentially want to use this way. Curious to hear thoughts from anyone on where it might come in handy.

[GuySten]

One case where it is potentialy useful to transport particles to nearest mesh boundary is when using weight windows and rouletting particles when crossing mesh surfaces.
I think I remember some discussion about that maybe in the discourse forum but I can't find it now...
Anyway according to the manual MCNP has an option to adjust particles weights when crossing weight windows boundaries.
I am not sure which cases that option is good for but we also might want it in openmc.

We can also potentially replace #3107 with adding a "timelike" mesh to stop particles at.
It might result in a cleaner code with less lines.

[JoffreyDorville]

Thanks for your comments @paulromano!

• I just updated the event-based mode as suggested.

• Regarding p.event_advance(), I like that all similar functions have the same signature (e.g., p.event_cross_surface(), p.event_collide()). But next_event is always used in the same context as p.event_advance() so this would save us a Particle attribute. I can change that if that's what you prefer here.

• I will probably focus on the temperature mesh for this PR, but I will keep the implementation as open as I can to allow for generalization later.

[paulromano]

A few initial comments; more to come!

[paulromano]

If the mesh used for the temperature field was used for something else and already written to XML, this would result in <values> never being written. Even if the mesh is shared, the values still need to be written.

[paulromano]

It would be good to have some validation here to check for cases that are not yet supported (unstructured mesh and, if I'm not mistaken, cylindrical/spherical mesh). Additionally, there could be some basic checks for consistency (number of values matches number of mesh elements, temperature values are positive).

[paulromano]

You should update docs/source/io_formats/settings.rst with any new XML elements.

[paulromano]

This should be declared in capi.h

[paulromano]

Mismatch between Doxygen comment and actual arguments (r vs bin)

[paulromano]

I see that ScalarField inherit from ABC but right now it doesn't have any abstract methods. Were you planning on making this class have a defined interface that subclasses would have to implement?

Also, right now these classes don't really have much "behavior" (methods), which calls into question if they are needed as classes in the first place, but I understand that they are likely to be expanded in the future. Can you elaborate on what you see these classes providing beyond what we could achieve with built-in type? i.e., why should it not just be:

settings.temperature_field = (mesh, values)