Fixed a bug in cell last update

src/particle.cpp

@@ -174,11 +174,6 @@ void Particle::event_calculate_xs()
   r_last() = r();
   time_last() = time();
 
-  // Reset event variables
-  event() = TallyEvent::KILL;
-  event_nuclide() = NUCLIDE_NONE;
-  event_mt() = REACTION_NONE;
-
   // If the cell hasn't been determined based on the particle's location,
   // initiate a search for the current cell. This generally happens at the
   // beginning of the history and again for any secondary particles
@@ -192,13 +187,18 @@ void Particle::event_calculate_xs()
     // Set birth cell attribute
     if (cell_born() == C_NONE)
       cell_born() = lowest_coord().cell();
+  }
 
-    // Initialize last cells from current cell
-    for (int j = 0; j < n_coord(); ++j) {
-      cell_last(j) = coord(j).cell();
-    }
-    n_coord_last() = n_coord();
+  // Initialize last cells from current cell
+  for (int j = 0; j < n_coord(); ++j) {
+    cell_last(j) = coord(j).cell();
   }
+  n_coord_last() = n_coord();
+
+  // Reset event variables
+  event() = TallyEvent::KILL;
+  event_nuclide() = NUCLIDE_NONE;
+  event_mt() = REACTION_NONE;
 
   // Write particle track.
   if (write_track())
@@ -232,6 +232,9 @@ void Particle::event_calculate_xs()
     macro_xs().fission = 0.0;
     macro_xs().nu_fission = 0.0;
   }
+
+  // Initialize last material from current material
+  material_last() = material();
 }
 
 void Particle::event_advance()
@@ -298,6 +301,9 @@ void Particle::event_cross_surface()
   }
   n_coord_last() = n_coord();
 
+  // Saving previous material data
+  material_last() = material();
+
   // Set surface that particle is on and adjust coordinate levels
   surface() = boundary().surface();
   n_coord() = boundary().coord_level();

tests/regression_tests/filter_cellfrom/results_true.dat

@@ -1,53 +1,53 @@
 k-combined:
 9.035025E-02 2.654309E-03
 tally 1:
-5.994069E+00
-3.594398E+00
+7.467961E+00
+5.580255E+00
 tally 2:
-4.559707E-04
-2.080739E-08
+0.000000E+00
+0.000000E+00
 tally 3:
-5.994525E+00
-3.594945E+00
+7.467961E+00
+5.580255E+00
 tally 4:
 0.000000E+00
 0.000000E+00
 tally 5:
-1.473892E+00
-2.187509E-01
+0.000000E+00
+0.000000E+00
 tally 6:
-5.223875E-05
-2.992861E-10
+7.821638E-04
+6.132129E-08
 tally 7:
-1.473945E+00
-2.187663E-01
+7.821638E-04
+6.132129E-08
 tally 8:
-1.885798E+01
-3.558423E+01
+0.000000E+00
+0.000000E+00
 tally 9:
 7.467961E+00
 5.580255E+00
 tally 10:
-5.082094E-04
-2.584432E-08
+7.821638E-04
+6.132129E-08
 tally 11:
-7.468470E+00
-5.581014E+00
+7.468743E+00
+5.581423E+00
 tally 12:
-1.885798E+01
-3.558423E+01
+0.000000E+00
+0.000000E+00
 tally 13:
 0.000000E+00
 0.000000E+00
 tally 14:
-2.739543E-04
-7.600983E-09
+0.000000E+00
+0.000000E+00
 tally 15:
-2.739543E-04
-7.600983E-09
+0.000000E+00
+0.000000E+00
 tally 16:
-7.881296E+01
-6.221087E+02
+9.767094E+01
+9.547828E+02
 tally 17:
 1.051397E+02
 1.106292E+03

tests/unit_tests/test_energy_balance.py

@@ -0,0 +1,105 @@
+import openmc
+import pytest
+
+
+@pytest.fixture
+def two_cell_model():
+    """Simple two-cell slab model with a fixed source.
+       Cell1 occupies x in [-10, 0], cell2 x in [0, 10]. 
+    """
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+
+    m1 = openmc.Material()
+    m1.add_element('Li', 1.0)
+    m1.set_density('g/cm3', 1.0)
+
+    m2 = openmc.Material()
+    m2.add_element('Al', 1.0)
+    m2.set_density('g/cm3', 1.0)
+
+    xmin = openmc.XPlane(-10.0, boundary_type="reflective")
+    xmid = openmc.XPlane(0.0)
+    xmax = openmc.XPlane(10.0, boundary_type="reflective")
+    ymin = openmc.YPlane(-10.0, boundary_type="reflective")
+    ymax = openmc.YPlane(10.0, boundary_type="reflective")
+    zmin = openmc.ZPlane(-10.0, boundary_type="reflective")
+    zmax = openmc.ZPlane(10.0, boundary_type="reflective")
+    cell1 = openmc.Cell(fill=m1, region=+xmin & -xmid & +ymin & -ymax & +zmin & -zmax)
+    cell2 = openmc.Cell(fill=m2, region=+xmid & -xmax & +ymin & -ymax & +zmin & -zmax)
+    model.geometry = openmc.Geometry([cell1, cell2])
+
+    src = openmc.IndependentSource()
+    src.space = openmc.stats.Point((-5.0, 0.0, 0.0))
+    src.particle = 'photon'
+    src.energy = openmc.stats.Discrete([1e6],[1.0])
+
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 1
+    model.settings.particles = 100
+    model.settings.source = src
+
+    return model, xmid, cell1, cell2, m1, m2
+
+
+def test_energy_balance_simple(two_cell_model, run_in_tmpdir):
+    model, xmid, cell1, cell2, *_  = two_cell_model
+
+    tally1 = openmc.Tally()
+    tally1.scores = ['heating']
+    tally1.filters = [openmc.CellFilter([cell1, cell2])]
+
+    tally2 = openmc.Tally()
+    tally2.scores = ['current']
+    tally2.filters = [openmc.EnergyFunctionFilter([0.0,20e6],[0.0,20e6]), openmc.SurfaceFilter([xmid])]
+
+
+    model.tallies = [tally1, tally2]
+    model.run(apply_tally_results=True)
+
+    assert tally1.mean.sum() == pytest.approx(1e6)
+
+    assert tally2.mean[0] == pytest.approx(tally1.mean[1])
+
+def test_current_conservation(two_cell_model, run_in_tmpdir):
+    model, xmid, cell1, cell2, m1, m2 = two_cell_model
+
+    energyfunc_filter = openmc.EnergyFunctionFilter([0.0,20e6],[0.0,20e6])
+
+    tally1 = openmc.Tally()
+    tally1.scores = ['current']
+    tally1.filters = [energyfunc_filter, openmc.SurfaceFilter([xmid])]
+
+    tally2 = openmc.Tally()
+    tally2.scores = ['current']
+    tally2.filters = [energyfunc_filter,
+                      openmc.CellFromFilter([cell1]),
+                      openmc.CellFilter([cell2])]
+
+    tally3 = openmc.Tally()
+    tally3.scores = ['current']
+    tally3.filters = [energyfunc_filter,
+                      openmc.CellFromFilter([cell2]),
+                      openmc.CellFilter([cell1])]
+
+    tally4 = openmc.Tally()
+    tally4.scores = ['current']
+    tally4.filters = [energyfunc_filter,
+                      openmc.MaterialFromFilter([m1]),
+                      openmc.MaterialFilter([m2])]
+
+    tally5 = openmc.Tally()
+    tally5.scores = ['current']
+    tally5.filters = [energyfunc_filter,
+                      openmc.MaterialFromFilter([m2]),
+                      openmc.MaterialFilter([m1])]
+
+
+    model.tallies = [tally1, tally2, tally3, tally4, tally5]
+    model.run(apply_tally_results=True)
+
+    assert pytest.approx(tally1.mean.sum()) == tally2.mean[0] + tally3.mean[0]
+
+    assert tally2.mean[0] == pytest.approx(tally4.mean[0])
+
+    assert tally3.mean[0] == pytest.approx(tally5.mean[0])