pshriwise commented 6 months ago

Short Version

A constant defining the structure of the EntityHandle space in MOAB was changed between MOAB 5.3.0 and MOAB 5.4.0 (commit linked below), revealing an issue with placement of the implicit complement handle in the DAGMC indices structure.

For a file with enough MeshSet's to be read, the implicit complement handle may not appear at the end of the EntityHandle space, meaning that one may perform a point containment check on the implicit complement volume before other, explicit volumes if looping over all volumes in the model. While its surprising that this constant changed, this is not the fault of MOAB, EntityHandles are required to be unique but not monotonically increasing when generated. This bug was present before the change in MOAB, but the change makes the likelihood of it occurring much higher.

Searching the implicit complement before other volumes can cause problems for source points that are coincident with surfaces. This may only be true for problems with imperfect topology relationships (failed imprinting/merging) and is why it didn't show up in the OpenMC tests, which use only well-formed models.

The fix for us would be to take extra care to ensure that the implicit complement volume is placed at the back of our data structures in DagMC::build_indices.

In external code interfaces, it would probably be wise to have a check that the implicit complement is checked last. In OpenMC, for example, this would look like the following in openmc::DAGUniverse::find_cell:

bool DAGUniverse::find_cell(Particle& p) const
{
  // if the particle isn't in any of the other DagMC
  // cells, place it in the implicit complement
  for (auto cell_idx : cells_) {
    // skip the implicit complement for now
    if (cell_idx == this->implicit_complement_idx()) continue;

    const auto& cell = model::cells[cell_idx];
    if(cell->contains(p.r(), p.u(), p.surface())) {
      p.lowest_coord().cell = cell_idx;
      return true;
    }
  }

  // finally, check the implicit complement for containment
  const auto& ipc_cell = model::cells[implicit_complement_idx()];
  if (cell->contains(p.r(), p.u(), p.surface())) {
    p.lowest_coord().cell = implicit_complement_idx();
    return true;
  }

  if (model::universe_map[this->id_] != model::root_universe) {
    p.lowest_coord().cell = implicit_complement_idx();
    return true;
  }
  return false;
}

Long (Long) Version

What changed in MOAB?

Between versions 5.3.0 and 5.4.0, a fundamental constant was changed in MOAB, SequenceManager::DEFAULT_VERTEX_SEQUENCE_SIZE , in this commit which propagates up to the variable SequenceManager::DEFAULT_MESHSET_SEQUENCE_SIZE. This update causes the sequence of EntityHandle's to change when generating new MeshSet's. This changes affects other entity types as well, but only the MeshSet EntityHandle's impact the problem here.

From debugging the GeomTopoTool::generate_implicit_complement method:

Under the call to create a new meshset for the implicit complement we get to SequenceManager::create_mesh_set. In both versions of MOAB we then call TypeSequenceManager::find_free_handle, presumably searching for a free entity handle in the space of the existing handle sequences that have already been allocated.
Inside TypeSequenceManager::find_free_handle, one loop goes over a std::set of

The difference between the two becomes evident when loading a file:

MOAB generates and pre-allocates memory for new handles in chunks, referred to as sequences, with a default size based on the variables changed in the commit linked above. The allocation is managed by a SequenceData object and associated EntitySequences indicate which EntityHandles have been assigned that sequence (via Core::create_meshset, Core::create_vertex, etc.). More on this topic can be found in the documentation here. When a new EntityHandle is requested and the current sequence is out of space, another is allocated automatically. New sequences are almost always generated with the default size, with some exceptions. File I/O is one of these exceptions and has an impact on loading files with DAGMC.

We create a new file set MeshSet in DagMC::load_file, this file gets the first EntityHandle available in the sequence (something like 12682136550675316737) in both MOAB versions. Nothing new. The behavior diverges when we load a file in ReadUtil::read_sets. This method looks for a MeshSetSequence that will hold all of the sets in the file. If one isn't found that will contain them all, then a new sequence is created behind the initial MeshSetSequence that was created when the file set was made. If a sequence exists that will hold all of the MeshSet's in the file, it is used starting with the next available EntityHandle in the sequence. In 5.3.x, the default EntityHandle sequence is holds 524,287 handles whereas in 5.4.1 the default size sequence holds 16,383 handles. The outcome is this:

5.3.1 EntityHandle Space: | File Set | .h5m file sets |

5.4.1 EntityHandle Space: | File Set | Unused Handle Space | .h5m file sets |

Later, when DAGMC creates the implicit complement meshset via moab::GeomTopoTool::generate_implicit_complement, the next available handle in the global sequence is used:

5.3.1 EntityHandle Space After Generating the Implicit Complement: | File Set | .h5m file sets | IPC Handle |

How does this affect DAGMC?

When the implicit complement is generated by DAGMC, a new MeshSet is created in the MOAB instance. Previously, the EntityHandle of this MeshSet was larger than all existing MeshSet's in the MOAB instance (at least in many/most cases??) and, as a result, this handle would end up at the end of the GeomTopoTool::geomRanges data member (a series of containers that will always iterate over the contained EntityHandle's in increasing order of their value) when iterating through the volume EntityHandle's. This ordering propagates into the DagMC::entHandles data member that is used by Monte Carlo code interfaces to create cells.

The effective change in behavior is that the implicit complement volume is checked for containment before other explicit volumes. In a model where all surfaces aren't merged or other problems exist, this can cause ambiguity in point containment queries. This was previously (perhaps intentionally?) resolved by favoring containment in explicitly defined volumes over the implicit complement volume. With the new ordering of EntityHandle's resulting from the upstream change in MOAB, this is no longer the case. So some particles are being born in regions of the implicit complement

Evidence of this problem in the Open MSRE Model

This is a model kindly generated and curated by Copenhagen Atomics here . @paulromano discovered a large discrepancy in the eigenvalue when doing nothing but changing between versions of MOAB in his software stack.

OpenMC k-eff w/ 5.3.0: 1.01285 +/- 0.00225 OpenMC k-eff w/ 5.4.1: 0.98537 +/- 0.00220

When debugging the problem and comparing tracks I noted that some tracks of identical position and direction were born in different cells. These cells had different IDs and different indexes in the DAGMC datastructures due to the new EntityHandle sequencing, making it difficult to determine the root cause of the issue. Eventually I discovered that the EntityHandle space of volumes were very different and that the implicit complement handle was at the end of the data structure rather than the front for MOAB 5.4.1:

MOAB 5.3.0

Pasted image 20240109134657

MOAB 5.4.1 Pasted image 20240109134718

pshriwise commented 5 months ago

Noting here that this affects other operations in DAGMC such as implicit complement material assignment, which assumes that we read the IPC's material assignment off a graveyard volume before getting to the IPC handle in this method.

https://github.com/svalinn/DAGMC/blob/b001729e89f9a8b6f958ea6bf3628319901c88a9/src/dagmc/dagmcmetadata.cpp#L112-L249

gonuke commented 5 months ago

Noting here that this affects other operations in DAGMC such as implicit complement material assignment, which assumes that we read the IPC's material assignment off a graveyard volume before getting to the IPC handle in this method.

https://github.com/svalinn/DAGMC/blob/b001729e89f9a8b6f958ea6bf3628319901c88a9/src/dagmc/dagmcmetadata.cpp#L112-L249

Does this mean we need a fix for this? Or does your fix take care of this?

pshriwise commented 5 months ago

935 should fix this as well. Just noting that our dependence on the IPC being at the back of the DAGMC index space has affects outside of the geometry queries.