Skip to content

Reevaluate cell in a close collision after a surface crossing#3933

Open
GuySten wants to merge 1 commit intoopenmc-dev:developfrom
GuySten:collision-cell-change
Open

Reevaluate cell in a close collision after a surface crossing#3933
GuySten wants to merge 1 commit intoopenmc-dev:developfrom
GuySten:collision-cell-change

Conversation

@GuySten
Copy link
Copy Markdown
Contributor

@GuySten GuySten commented May 3, 2026

Description

This PR adds a cell lookup after a collision that could change particle's cell.
This PR is an alternative path for #3923.

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 GuySten requested a review from paulromano May 3, 2026 08:06
@GuySten
Copy link
Copy Markdown
Contributor Author

GuySten commented May 3, 2026

@lewisgross1296, can you check if this PR improves your problem as well?

@GuySten GuySten changed the title recalc cell in a collision after a surface crossing Reevaluate cell in a close collision after a surface crossing May 3, 2026
@GuySten GuySten added the Bugs label May 3, 2026
@lewisgross1296
Copy link
Copy Markdown
Contributor

lewisgross1296 commented May 3, 2026

@lewisgross1296, can you check if this PR improves your problem as well?

Just rebuilt with this PR and launched my openmc -g. Probably will take at least one day to know for sure. Thanks for the quick alternative!

@lewisgross1296
Copy link
Copy Markdown
Contributor

Checked back this morning. Running the same XML with the code in this branch (and some extra prints in src/geometry.cpp) led to the same failure mode I saw before. The output below shows the particle detects two overlaps before encountering the maximum number of events and failing to bank a fission site

       56/1    1.01606    1.01589 +/- 0.00012
Logical checks
currently on coordinate level 2
p.coord(j).r(), p.coord(j).u() = (-0.456426, 0.553288, -1.60396),(-0.79841, 0.0776792, 0.635684)
index_cell =566
p.coord(j).cell() 565
cell instance 56
Surface token for the surface that the particle is currently on 0
GeometryState information
r_born = (9.15817, 211.291, 10.351)
Current global coordinate position r=(18.6115, 228.396, 70.5)
Current global coordinate direction u=(-0.79841, 0.635684, -0.0776792)
Coordinates of last collision or reflective/periodic surface crossing for current tallies(18.6115, 228.396, 70.5)
Previous spatial coordinates before collision r=(18.6115, 228.396, 70.5)
Previous direction vector before last collision u=(-0.79841, 0.635684, -0.0776792)
Current material 1
Last material -1
printing cell info
distribcell_index_234
Overlapping cells detected: 671, 670 on universe 346

Logical checks
currently on coordinate level 2
p.coord(j).r(), p.coord(j).u() = (-0.731418, 0.580042, -1.38501),(-0.310974, 0.523344, 0.814965)
index_cell =566
p.coord(j).cell() 565
cell instance 56
Surface token for the surface that the particle is currently on 0
GeometryState information
r_born = (9.15817, 211.291, 10.351)
Current global coordinate position r=(18.3365, 228.615, 70.4733)
Current global coordinate direction u=(-0.310974, 0.814965, -0.523344)
Coordinates of last collision or reflective/periodic surface crossing for current tallies(18.3365, 228.615, 70.4733)
Previous spatial coordinates before collision r=(18.3365, 228.615, 70.4733)
Previous direction vector before last collision u=(-0.310974, 0.814965, -0.523344)
Current material 1
Last material -1
printing cell info
distribcell_index_234
Overlapping cells detected: 671, 670 on universe 346
 WARNING: Particle 2227797 underwent maximum number of events.
       57/1   -nan       -nan     +/- 0.00000
 ERROR: No fission sites banked on MPI rank 0

The trouble is happening on the same batch as before and is "detecting" an overlap on the same two cells. Regarding @paulromano's diagnosis in #3929, the above might be related to the case he described.

In the printout, Current material 1 is my helium coolant (density 0.00295267 g/cm^3). However, the particle overlap is detected in material 2 which is my reflector (density 2.7987 g/cm^3). The overlapping cells (670,671) are both made of Material 2. The current global coordinate position is in cell 671 (material 2), however, if I understand the logic correctly, the actual cell it's in is 670 (model::cells[p.coord(j).cell()]->id_ = 670).

This location is near a coolant channel, but I'm a little confused why it thinks it's in Current material 1, as the global coordinate position is in material 2. As discussed with @gonuke in the past, we think some event may have caused the particle's actual position to become displaced from the position it thinks it's in.

I uploaded my XML on discourse. The printouts above in come from modifying geometry.cpp

bool check_cell_overlap(GeometryState& p, bool error)
{
  int n_coord = p.n_coord();
  // Loop through each coordinate level
  for (int j = 0; j < n_coord; j++) {
    Universe& univ = *model::universes[p.coord(j).universe()];

    // Loop through each cell on this level
    for (auto index_cell : univ.cells_) {
      Cell& c = *model::cells[index_cell];
      if (c.contains(p.coord(j).r(), p.coord(j).u(), p.surface())) {
        if (index_cell != p.coord(j).cell()) {
          if (error) {
            // Print info of check
            std::cout << "Logical checks" << std::endl;
            std::cout << "currently on coordinate level " << j << std::endl;
            std::cout << "p.coord(j).r(), p.coord(j).u() = " << p.coord(j).r() << 
            "," << p.coord(j).u() << std::endl;
            std::cout << "index_cell =" << index_cell << std::endl;
            std::cout <<  "p.coord(j).cell() " << p.coord(j).cell() << std::endl;
            std::cout << "cell instance " << p.cell_instance() << std::endl;
            std::cout << "Surface token for the surface that the particle is currently on " << p.surface() << std::endl;
            std::cout << "" << std::endl;

            // PRINT INFO
            std::cout << "GeometryState information" << std::endl;
            std::cout << "r_born = " << p.r_born() << std::endl;
            std::cout << "Current global coordinate position r=" << p.r() << std::endl;
            std::cout << "Current global coordinate direction u=" << p.u() << std::endl;
            std::cout << "Coordinates of last collision or reflective/periodic surface crossing for current tallies" << p.r_last_current() << std::endl;
            std::cout << "Previous spatial coordinates before collision r=" << p.r_last() << std::endl;
            std::cout << "Previous direction vector before last collision u=" << p.u_last() << std::endl;
            std::cout << "Current material " << p.material() << std::endl;
            std::cout << "Last material " << p.material_last() << std::endl;

            // cell class info
            std::cout << "printing cell info" << std::endl;
            std::cout << "distribcell_index_" << c.distribcell_index_ << std::endl;

            auto s = fmt::format("Overlapping cells detected: {}, {} on universe {}",
                c.id_, model::cells[p.coord(j).cell()]->id_, univ.id_);
            std::cout << s << std::endl;
            std::cout << std::endl;

            // turn off fatal error and see if we keep getting overlaps
            // fatal_error(
            //   fmt::format("Overlapping cells detected: {}, {} on universe {}",
            //     c.id_, model::cells[p.coord(j).cell()]->id_, univ.id_));
          }
          return true;
        }
#pragma omp atomic
        ++model::overlap_check_count[index_cell];
      }
    }
  }

  return false;
}

@GuySten
Copy link
Copy Markdown
Contributor Author

GuySten commented May 4, 2026

@lewisgross1296, I think your problem may come from something other that small cylinders.
The problem is that there are a lot of potential candidates for what is going on:

  1. Maybe there is a problem with universe transformations (your example case has complex nesting).
  2. Some event caused the particle position to be displaced.
  3. The particle is close to a surface boundary and it's direction is perpendicular to the surface normal (So it is not clear if it is entering or exiting the cell)
  4. Something else...
  5. A combination of the above.
    Digging into this problem is very hard without a minimal example (that is reproducible and completes in a reasonable amout of time).
    I think for the time being a practical solution will be to change your seed a few times until your problem succeeds.

@lewisgross1296
Copy link
Copy Markdown
Contributor

Digging into this problem is very hard without a minimal example (that is reproducible and completes in a reasonable amout of time).

I agree. This has been really tough to make progress on due to how long it takes it test anything / how rarely this occurs. Unfortunately, I need to use a design this complicated. If I could make a more minimal example, that would really help, but I'm not totally sure I actually understand the root issue in a way that would allow me to confidently create an MWE.

I have heard that the intersection of planes and other surface types (sphere, cylinder) can lead to particles getting displaced from their true position (@gonuke called this "tunneling"). We do have many cylinders sliced by the planes that correspond to 1/6th symmetry, so perhaps that is leading to the problem. Strangely, I do not have the same issue with reflective boundary conditions (no overlaps at all despite the same high particle count), so my plan is to use those for now until we can try to figure out what is really going on here.

@GuySten GuySten requested review from pshriwise and shimwell May 5, 2026 09:49
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants