Determine surface-source starting half-space from particle direction

[paulromano]

Description

Recently, I've been working with some colleagues on a two-step workflow for a shielding calculation. In the first step, particles are accumulated on a surface through surface source, and in the second step the surface source is used as the starting source. We ran into problems with particles getting lost in the second step and I've traced it down to some bugs with how surface IDs that are written to a surface source file are used when that source file is used as the source in a new calculation. First, the surface ID that is written to the source file is unsigned, but when OpenMC reads the source file it treats the surface ID as being signed (i.e., always on the positive side of the surface). Second, if the surface ID is not present in the model, rather than being ignored it ends up defaulting to the first surface in the model by virtue of calling surface_map[surf_id] when surf_id is not in the map.

Rather than write a signed surface ID to the surface source file, I've fixed the underlying issue by having surface source particles determine their starting half-space using the particle direction and the referenced surface normal. The surface ID is used only when the ID exists and the source position is coincident with that CSG surface; missing or mismatched IDs fall back to normal cell search.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) 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)

[GuySten]

Except for one suggestion.
Looks good to me.

[GuySten]

Suggested change

	if (site.surf_id != SURFACE_NONE) {
	auto it = model::surface_map.find(std::abs(site.surf_id));
	if (it == model::surface_map.end()) {
	site.surf_id = SURFACE_NONE;
	} else {
	const auto& surf = *model::surfaces[it->second];
	if (surf.geom_type() == GeometryType::CSG &&
	std::abs(surf.evaluate(site.r)) < FP_COINCIDENT) {
	int surf_id = std::abs(site.surf_id);
	site.surf_id =
	(site.u.dot(surf.normal(site.r)) > 0.0) ? surf_id : -surf_id;
	} else {
	site.surf_id = SURFACE_NONE;
	}
	}
	}
	return site;
	if (site.surf_id != SURFACE_NONE) {
	auto it = model::surface_map.find(std::abs(site.surf_id));
	if (it != model::surface_map.end()) {
	const auto& surf = *model::surfaces[it->second];
	if (surf.geom_type() == GeometryType::CSG &&
	std::abs(surf.evaluate(site.r)) < FP_COINCIDENT) {
	int surf_id = std::abs(site.surf_id);
	site.surf_id =
	(site.u.dot(surf.normal(site.r)) > 0.0) ? surf_id : -surf_id;
	return site;
	}
	}
	site.surf_id = SURFACE_NONE;
	}
	return site;