Abstract the geometry module

[smharper]

I want to add a layer of abstraction on top of the geometry module so we can easily implement different geometry models. For example, @friedmud has expressed interest in using finite element geometry directly in OpenMC to make coupling with multiphysics simulations easier. I am also interested in tracking on boundary representation (b-rep) geometry simply because we could use a lot of existing software to easily model complex geometries (I mostly have Blender in mind here).

An abstract base class would mean that a developer who is working on an alternate geometry model won't have to muck about with global memory like cells and the rather overloaded particle object. They might even be able to glue on existing geometry libraries that perhaps aren't written in Fortran.

Here's what I'm thinking for the abstraction:

We'll have an abstract base class called GeometryKernel. All of the data related to the location of particles will be in memory controlled exclusively by the kernel. In other words, the linked list of coords will no longer be part of the particle data structure; it will be stored somewhere in the geometry kernel. Also, all of the geometry stuff like the arrays of cells, universes, etc. will be moved into the kernel. It will have to provide some bookkeeping methods like this:

• subroutine initialize(): Reads the XML file, initializes internal arrays, stuff like that.
• subroutine allocate_particle(integer particle_id, real*8 xyz(3), real*8 uvw(3)): Prepares the kernel memory for one particle (like making a coord array). The rest of OpenMC will only be able to access geometry data about that particle by requesting it from the kernel, using the particle_id number.
• subroutine deallocate_particle(integer particle_id)

In order to store source sites and handle scattering reactions we'll need:

• function get_xyz(integer particle_id) result(real*8 xyz(3))
• function get_uvw(integer particle_id) result(real*8 uvw(3))
• subroutine set_uvw(integer particle_id, real*8 uvw(3))

For tracking:

• function get_distance_to_boundary(integer particle_id) result(real*8 distance)
• subroutine translate_to_boundary(integer particle_id)

Simple tallying stuff:

• function get_cell_ids(integer particle_id result(integer cell_ids(:)): This will return an array of ids for all of the cells that the particle is in.
• function get_unique_cell_ids(integer particle_id) result(integer cell_uids(:)): This would replace the offset maps we use now for distribcell.
• function get_surface_id(integer particle_id) result(integer surf_id): For surface current tallies. I imagine it will only be valid if translate_to_boundary has just been called.

Mesh tallies: After some internal debate, I've decided that all of the guts of the mesh tallies should remain in the tally module. This means that the tracking module will have to remember the initial and final track coordinates from get_xyz and pass them into the tally module.

Functional expansion tallies: We don't have these yet (at least, not on the xyz basis), but they are likely coming so we need to plan for them:

• function get_local_xyz(integer particle_id, integer unique_cell_id) result(real*8 xyz(3))

I'm not decided on how to handle errors. We could add an error argument like integer, intent(inout) :: i_err to all of the methods, but I think that's really ugly. Any suggestions?

One last issue I'd like to mention is how to handle materials. For the near-term, I think it would be fine for the tracking module to have a map from cell ids to materials. But eventually, I think the geometry module should be providing the material information directly. This will make it easier for us to handle distribmats (which we are working on currently, right?) and continuously varying materials (which at least some of us are dreaming about). But I have no idea how to handle that so I'll leave it for another time.

[friedmud]

:+1:

[paulromano]

Some abstraction is badly needed on the geometry side, so I'm glad you've been thinking about this @smharper. Aside from having geometry kernels that follow the same type of interface, are you planning on abstracting the CSG geometry itself as part of this? For example, having a Surface base class that is inherited by different surface types.

My one comment is regarding moving the particle coordinates to the geometry module. The current setup is very 'clean' for multithreading because each particle is already threadprivate. Moving it elsewhere means that the Particle class no longer contains all the threadprivate information related to the particle. What is your reasoning for wanting to move the coordinates to geometry, rather than having the API calls include an argument for the coordinates?

[smharper]

I'm not planning on abstracting the CSG primitives right now, but I agree that it's a good idea, and I would like to get it done eventually.

I want to move the coordinates to the geometry kernel for two reasons: First, I think they are too specific to CSG. BRep or FE geometry will likely not have this concept of universes, nested levels, and lattices. However, those geometry models would likely need k-d tree information which would also need to be threadprivate. Second, it makes the code modules less tightly coupled together, which in turn makes it easier to change the geometry implementation without messing up other sections of the code. A developer who is modifying the geometry kernel wouldn't have to worry about how their changes to the implementation might cause errors in the tally, tracking, plotting, etc. modules.

[paulromano]

All very good points. Thanks for clarifying!

[gridley]

I would argue that @pshriwise's abstraction of the Universe base class resolves this issue. BRep could be its own universe or FE meshes could also be a universe, which DAGUniverse shows how to do.

The threadprivateness I do not think is relevant to consider in closing this. This acceleration data would likely be implemented as a custom partitioner as illustrated by Saad's PR #2593, and if it requires extra geometry state, this would be tacked onto the GeometryState class as RayHistory is.