Trying to fix cell last

src/geometry.cpp

@@ -173,12 +173,21 @@ bool find_cell_inner(
       }
 
       // Set the material, temperature and density multiplier.
-      p.material_last() = p.material();
-      p.material() = c.material(p.cell_instance());
-      p.sqrtkT_last() = p.sqrtkT();
-      p.sqrtkT() = c.sqrtkT(p.cell_instance());
-      p.density_mult_last() = p.density_mult();
-      p.density_mult() = c.density_mult(p.cell_instance());
+      if (p.sqrtkT() != -1) {
+        // if current material is set, set last material
+        p.material_last() = p.material();
+        p.material() = c.material(p.cell_instance());
+        p.sqrtkT_last() = p.sqrtkT();
+        p.sqrtkT() = c.sqrtkT(p.cell_instance());
+        p.density_mult_last() = p.density_mult();
+        p.density_mult() = c.density_mult(p.cell_instance());
+      } else {
+        // set both current data and last data to current data
+        p.material_last() = p.material() = c.material(p.cell_instance());
+        p.sqrtkT_last() = p.sqrtkT() = c.sqrtkT(p.cell_instance());
+        p.density_mult_last() = p.density_mult() =
+          c.density_mult(p.cell_instance());
+      }
 
       return true;
 

src/particle.cpp

@@ -198,6 +198,8 @@ void Particle::event_calculate_xs()
       cell_last(j) = coord(j).cell();
     }
     n_coord_last() = n_coord();
+
+    material_last() = material();
   }
 
   // Write particle track.
@@ -211,7 +213,7 @@ void Particle::event_calculate_xs()
   if (material() != MATERIAL_VOID) {
     if (settings::run_CE) {
       if (material() != material_last() || sqrtkT() != sqrtkT_last() ||
-          density_mult() != density_mult_last()) {
+          density_mult() != density_mult_last() || true) {
         // If the material is the same as the last material and the
         // temperature hasn't changed, we don't need to lookup cross
         // sections again.
@@ -397,10 +399,6 @@ void Particle::event_collide()
   // Save coordinates for tallying purposes
   r_last_current() = r();
 
-  // Set last material to none since cross sections will need to be
-  // re-evaluated
-  material_last() = C_NONE;
-
   // Set all directions to base level -- right now, after a collision, only
   // the base level directions are changed
   for (int j = 0; j < n_coord() - 1; ++j) {
@@ -419,6 +417,14 @@ void Particle::event_collide()
   if (!model::active_tallies.empty())
     score_collision_derivative(*this);
 
+  // Saving previous cell data
+  for (int j = 0; j < n_coord(); ++j) {
+    cell_last(j) = coord(j).cell();
+  }
+  n_coord_last() = n_coord();
+
+  material_last() = material();
+
 #ifdef OPENMC_DAGMC_ENABLED
   history().reset();
 #endif
@@ -472,6 +478,8 @@ void Particle::event_revive_from_secondary()
           cell_last(j) = coord(j).cell();
         }
         n_coord_last() = n_coord();
+
+        material_last() = material();
       }
       pht_secondary_particles();
     }
@@ -696,6 +704,7 @@ void Particle::cross_reflective_bc(const Surface& surf, Direction new_u)
 
   // Reassign particle's cell and surface
   coord(0).cell() = cell_last(0);
+  material() = material_last();
   surface() = -surface();
 
   // If a reflective surface is coincident with a lattice or universe

tests/regression_tests/filter_cellfrom/results_true.dat

@@ -1,29 +1,29 @@
 k-combined:
 9.035025E-02 2.654309E-03
 tally 1:
-5.994069E+00
-3.594398E+00
+6.586649E+00
+4.340013E+00
 tally 2:
 4.559707E-04
 2.080739E-08
 tally 3:
-5.994525E+00
-3.594945E+00
+6.587105E+00
+4.340613E+00
 tally 4:
 0.000000E+00
 0.000000E+00
 tally 5:
-1.473892E+00
-2.187509E-01
+8.813120E-01
+7.829296E-02
 tally 6:
 5.223875E-05
 2.992861E-10
 tally 7:
-1.473945E+00
-2.187663E-01
+8.813643E-01
+7.830214E-02
 tally 8:
-1.885798E+01
-3.558423E+01
+5.908648E+00
+3.492172E+00
 tally 9:
 7.467961E+00
 5.580255E+00
@@ -34,8 +34,8 @@ tally 11:
 7.468470E+00
 5.581014E+00
 tally 12:
-1.885798E+01
-3.558423E+01
+5.908648E+00
+3.492172E+00
 tally 13:
 0.000000E+00
 0.000000E+00
@@ -46,8 +46,8 @@ tally 15:
 2.739543E-04
 7.600983E-09
 tally 16:
-7.881296E+01
-6.221087E+02
+9.176229E+01
+8.428689E+02
 tally 17:
 1.051397E+02
 1.106292E+03

tests/unit_tests/test_energy_balance.py

@@ -0,0 +1,152 @@
+import openmc
+import pytest
+import numpy as np
+
+
+@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
+
+@pytest.fixture
+def two_sphere_model():
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+
+    mat = openmc.Material()
+    mat.add_nuclide('H2', 1.0)
+    mat.set_density('g/cm3', 1.0)
+
+    surf1 = openmc.Sphere(r=0.01)
+    surf2 = openmc.Sphere(r=1000, boundary_type='vacuum')
+    cell1 = openmc.Cell(region=-surf1)
+    cell2 = openmc.Cell(fill=mat, region=+surf1 & -surf2)
+
+    model.geometry = openmc.Geometry([cell1, cell2])
+
+    src = openmc.IndependentSource()
+    src.energy = openmc.stats.Discrete([1e6],[1.0])
+
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 1
+    model.settings.particles = 1000
+    model.settings.source = src
+
+    return model, cell1, cell2, mat
+
+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])
+
+def test_cellfrom_heating(run_in_tmpdir, two_sphere_model):
+
+    model, cell1, cell2, mat = two_sphere_model
+
+    tally1 = openmc.Tally()
+    tally1.scores = ['heating']
+    tally1.filters = [openmc.CellFromFilter([cell1, cell2]), openmc.CellFilter([cell2])]
+
+    tally2 = openmc.Tally()
+    tally2.scores = ['heating']
+    tally2.filters = [openmc.MaterialFromFilter([mat]), openmc.MaterialFilter([mat])]
+
+    model.tallies = [tally1, tally2]
+
+    model.run(apply_tally_results=True)
+
+    assert np.all(tally1.mean > 0)
+
+    assert tally1.mean[1] == pytest.approx(tally2.mean[0]) 
+