Implicit Complement Position in DAGMC Indices

[pshriwise]

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 |

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

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

MOAB 5.4.1

[pshriwise]

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.

DAGMC/src/dagmc/dagmcmetadata.cpp

Lines 112 to 249 in b001729

	std::string implicit_complement_material{};
	std::string implicit_complement_density{};
	// loop over all cells
	for (int i = 1; i <= num_cells; ++i) {
	int cellid = DAG->id_by_index(3, i);
	moab::EntityHandle eh = DAG->entity_by_index(3, i);
	material_props = material_assignments[eh];
	density_props = density_assignments[eh];
	// this is actually ok for a single volume, one that has the _comp tag at
	// the end of it
	if (material_props.size() > 1) {
	// search the props for _comp
	std::size_t comp_found;
	int position;
	for (int j = 0; j < material_props.size(); j++) {
	comp_found = material_props[j].find("_comp");
	if (comp_found != std::string::npos) {
	position = j;
	break;
	}
	}
	if (comp_found != std::string::npos) {
	// success found the _comp tag for the impl_compl material
	// set the impl_comp material for use later
	implicit_complement_material = material_props[position].substr(
	0, material_props[position].size() - 5);
	implicit_complement_density = density_props[0];
	material_props.erase(material_props.begin() + position);
	} else {
	// failure, a volume can only have a single material associated with it
	std::stringstream ss;
	ss << "More than one material for volume with id " << cellid
	<< std::endl;
	ss << "that does not the the _comp tag associated with it."
	<< std::endl;
	ss << cellid << " has the following material assignments:" << std::endl;
	for (int j = 0; j < material_props.size(); j++) {
	ss << material_props[j] << std::endl;
	}
	ss << "Please check your material assignments " << cellid;
	logger.error(ss.str());
	exit(EXIT_FAILURE);
	}
	} else {
	// because of how the data are inserted into the map, there is always
	// at least one entry, "" if nothing is found
	// if there is no material property - not failure for impl_comp
	if (material_props[0] == "" && !(DAG->is_implicit_complement(eh))) {
	std::stringstream ss;
	ss << "No material property found for volume with ID " << cellid
	<< std::endl;
	ss << "Every volume must have exactly one mat: property";
	logger.error(ss.str());
	exit(EXIT_FAILURE);
	}
	}
	// this is never ok for a volume to have more than one property for density
	if (density_props.size() > 1) {
	std::stringstream ss;
	ss << "More than one density specified for " << cellid << std::endl;
	ss << cellid << " has the following density assignments" << std::endl;
	for (int j = 0; j < density_props.size(); j++) {
	ss << density_props[j] << std::endl;
	}
	ss << "Please check your density assignments " << cellid;
	logger.error(ss.str());
	exit(EXIT_FAILURE);
	}
	std::string grp_name{};
	// determine if we have a density property
	if (density_props.size() == 1) {
	if (!density_props[0].empty()) {
	grp_name = "mat:" + material_props[0] + "/rho:" + density_props[0];
	volume_density_data_eh[eh] = density_props[0];
	} else {
	grp_name = "mat:" + material_props[0];
	volume_density_data_eh[eh] = "";
	}
	}
	// check to see if the simplified naming scheme is used, by try to convert
	// the material property to an int
	if (require_density && try_to_make_int(material_props[0]) &&
	density_props[0].empty() && !(DAG->is_implicit_complement(eh))) {
	std::stringstream ss;
	ss << "Using the simplified naming scheme without a density" << std::endl;
	ss << "property is forbidden, please rename the group mat:"
	<< material_props[0];
	logger.error(ss.str());
	exit(EXIT_FAILURE);
	}
	// set the material value
	volume_material_property_data_eh[eh] = grp_name;
	bool is_graveyard =
	to_lower(grp_name).find(to_lower(graveyard_str)) != std::string::npos;
	bool is_vacuum =
	to_lower(grp_name).find(to_lower(vacuum_str)) != std::string::npos;
	// not graveyard or vacuum or implicit compliment
	if (!is_graveyard && !is_vacuum && !DAG->is_implicit_complement(eh)) {
	volume_material_data_eh[eh] = material_props[0];
	}
	// found graveyard
	else if (is_graveyard) {
	volume_material_property_data_eh[eh] = "mat:Graveyard";
	volume_material_data_eh[eh] = graveyard_str;
	}
	// vacuum
	else if (is_vacuum) {
	volume_material_property_data_eh[eh] = "mat:Vacuum";
	volume_material_data_eh[eh] = vacuum_str;
	}
	// implicit complement
	else if (DAG->is_implicit_complement(eh)) {
	if (implicit_complement_material == "") {
	logger.message("Implicit Complement assumed to be Vacuum");
	volume_material_property_data_eh[eh] = "mat:Vacuum";
	volume_material_data_eh[eh] = vacuum_str;
	} else {
	volume_material_property_data_eh[eh] =
	"mat:" + implicit_complement_material;
	if (implicit_complement_density != "")
	volume_material_property_data_eh[eh] +=
	"/rho" + implicit_complement_density;
	volume_material_data_eh[eh] = implicit_complement_material;
	volume_density_data_eh[eh] = implicit_complement_density;
	}
	}
	}
	}

[gonuke]

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.

DAGMC/src/dagmc/dagmcmetadata.cpp

Lines 112 to 249 in b001729

	std::string implicit_complement_material{};
	std::string implicit_complement_density{};
	// loop over all cells
	for (int i = 1; i <= num_cells; ++i) {
	int cellid = DAG->id_by_index(3, i);
	moab::EntityHandle eh = DAG->entity_by_index(3, i);
	material_props = material_assignments[eh];
	density_props = density_assignments[eh];
	// this is actually ok for a single volume, one that has the _comp tag at
	// the end of it
	if (material_props.size() > 1) {
	// search the props for _comp
	std::size_t comp_found;
	int position;
	for (int j = 0; j < material_props.size(); j++) {
	comp_found = material_props[j].find("_comp");
	if (comp_found != std::string::npos) {
	position = j;
	break;
	}
	}
	if (comp_found != std::string::npos) {
	// success found the _comp tag for the impl_compl material
	// set the impl_comp material for use later
	implicit_complement_material = material_props[position].substr(
	0, material_props[position].size() - 5);
	implicit_complement_density = density_props[0];
	material_props.erase(material_props.begin() + position);
	} else {
	// failure, a volume can only have a single material associated with it
	std::stringstream ss;
	ss << "More than one material for volume with id " << cellid
	<< std::endl;
	ss << "that does not the the _comp tag associated with it."
	<< std::endl;
	ss << cellid << " has the following material assignments:" << std::endl;
	for (int j = 0; j < material_props.size(); j++) {
	ss << material_props[j] << std::endl;
	}
	ss << "Please check your material assignments " << cellid;
	logger.error(ss.str());
	exit(EXIT_FAILURE);
	}
	} else {
	// because of how the data are inserted into the map, there is always
	// at least one entry, "" if nothing is found
	// if there is no material property - not failure for impl_comp
	if (material_props[0] == "" && !(DAG->is_implicit_complement(eh))) {
	std::stringstream ss;
	ss << "No material property found for volume with ID " << cellid
	<< std::endl;
	ss << "Every volume must have exactly one mat: property";
	logger.error(ss.str());
	exit(EXIT_FAILURE);
	}
	}
	// this is never ok for a volume to have more than one property for density
	if (density_props.size() > 1) {
	std::stringstream ss;
	ss << "More than one density specified for " << cellid << std::endl;
	ss << cellid << " has the following density assignments" << std::endl;
	for (int j = 0; j < density_props.size(); j++) {
	ss << density_props[j] << std::endl;
	}
	ss << "Please check your density assignments " << cellid;
	logger.error(ss.str());
	exit(EXIT_FAILURE);
	}
	std::string grp_name{};
	// determine if we have a density property
	if (density_props.size() == 1) {
	if (!density_props[0].empty()) {
	grp_name = "mat:" + material_props[0] + "/rho:" + density_props[0];
	volume_density_data_eh[eh] = density_props[0];
	} else {
	grp_name = "mat:" + material_props[0];
	volume_density_data_eh[eh] = "";
	}
	}
	// check to see if the simplified naming scheme is used, by try to convert
	// the material property to an int
	if (require_density && try_to_make_int(material_props[0]) &&
	density_props[0].empty() && !(DAG->is_implicit_complement(eh))) {
	std::stringstream ss;
	ss << "Using the simplified naming scheme without a density" << std::endl;
	ss << "property is forbidden, please rename the group mat:"
	<< material_props[0];
	logger.error(ss.str());
	exit(EXIT_FAILURE);
	}
	// set the material value
	volume_material_property_data_eh[eh] = grp_name;
	bool is_graveyard =
	to_lower(grp_name).find(to_lower(graveyard_str)) != std::string::npos;
	bool is_vacuum =
	to_lower(grp_name).find(to_lower(vacuum_str)) != std::string::npos;
	// not graveyard or vacuum or implicit compliment
	if (!is_graveyard && !is_vacuum && !DAG->is_implicit_complement(eh)) {
	volume_material_data_eh[eh] = material_props[0];
	}
	// found graveyard
	else if (is_graveyard) {
	volume_material_property_data_eh[eh] = "mat:Graveyard";
	volume_material_data_eh[eh] = graveyard_str;
	}
	// vacuum
	else if (is_vacuum) {
	volume_material_property_data_eh[eh] = "mat:Vacuum";
	volume_material_data_eh[eh] = vacuum_str;
	}
	// implicit complement
	else if (DAG->is_implicit_complement(eh)) {
	if (implicit_complement_material == "") {
	logger.message("Implicit Complement assumed to be Vacuum");
	volume_material_property_data_eh[eh] = "mat:Vacuum";
	volume_material_data_eh[eh] = vacuum_str;
	} else {
	volume_material_property_data_eh[eh] =
	"mat:" + implicit_complement_material;
	if (implicit_complement_density != "")
	volume_material_property_data_eh[eh] +=
	"/rho" + implicit_complement_density;
	volume_material_data_eh[eh] = implicit_complement_material;
	volume_density_data_eh[eh] = implicit_complement_density;
	}
	}
	}
	}

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

[pshriwise]

#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.

[gonuke]

Closed by #935