Add transformation boundary condition

include/openmc/boundary_condition.h

@@ -152,5 +152,24 @@ class RotationalPeriodicBC : public PeriodicBC {
   double angle_;
 };
 
+//==============================================================================
+//! A BC that transforms particle direction and position by matrices.
+//==============================================================================
+
+class TransformationBC : public BoundaryCondition {
+public:
+  TransformationBC(vector<double> dir_trans, vector<double> pos_trans);
+
+  void handle_particle(Particle& p, const Surface& surf) const override;
+
+  std::string type() const override { return "transformation"; }
+
+protected:
+  //! Matrix by which particle directions are transformed
+  vector<double> dir_trans_;
+  //! Matrix by which particle positions are transformed
+  vector<double> pos_trans_;
+};
+
 } // namespace openmc
 #endif // OPENMC_BOUNDARY_CONDITION_H

include/openmc/particle.h

@@ -103,6 +103,15 @@ class Particle : public ParticleData {
   void cross_periodic_bc(
     const Surface& surf, Position new_r, Direction new_u, int new_surface);
 
+  //! Cross a transformation boundary condition.
+  //
+  //! \param surf The surface (with the transformation boundary condition) that
+  //! the particle struck.
+  //! \param new_r The position of the particle after transformation.
+  //! \param new_u The direction of the particle after transformation.
+  void cross_transformation_bc(
+    const Surface& surf, Position new_r, Direction new_u);
+
   //! mark a particle as lost and create a particle restart file
   //! \param message A warning message to display
   virtual void mark_as_lost(const char* message) override;

src/boundary_condition.cpp

@@ -290,4 +290,41 @@ void RotationalPeriodicBC::handle_particle(
   p.cross_periodic_bc(surf, new_r, new_u, new_surface);
 }
 
+//==============================================================================
+// TransformationBC implementation
+//==============================================================================
+
+TransformationBC::TransformationBC(
+  vector<double> dir_trans, vector<double> pos_trans)
+{
+  // Set the transformation matrices
+  dir_trans_ = dir_trans;
+  pos_trans_ = pos_trans;
+}
+
+void TransformationBC::handle_particle(Particle& p, const Surface& surf) const
+{
+  Position r = p.r();
+  Direction u = p.u();
+
+  Position new_r = {r.x * pos_trans_[0] + r.y * pos_trans_[1] +
+                      r.z * pos_trans_[2] + pos_trans_[3],
+    r.x * pos_trans_[4] + r.y * pos_trans_[5] + r.z * pos_trans_[6] +
+      pos_trans_[7],
+    r.x * pos_trans_[8] + r.y * pos_trans_[9] + r.z * pos_trans_[10] +
+      pos_trans_[11]};
+
+  Position new_u = {u.x * dir_trans_[0] + u.y * dir_trans_[1] +
+                      u.z * dir_trans_[2] + dir_trans_[3],
+    u.x * dir_trans_[4] + u.y * dir_trans_[5] + u.z * dir_trans_[6] +
+      dir_trans_[7],
+    u.x * dir_trans_[8] + u.y * dir_trans_[9] + u.z * dir_trans_[10] +
+      dir_trans_[11]};
+
+  // Handle the effects of the surface albedo on the particle's weight.
+  BoundaryCondition::handle_albedo(p, surf);
+
+  p.cross_transformation_bc(surf, new_r, new_u);
+}
+
 } // namespace openmc

src/particle.cpp

@@ -750,6 +750,62 @@ void Particle::cross_periodic_bc(
   }
 }
 
+void Particle::cross_transformation_bc(
+  const Surface& surf, Position new_r, Direction new_u)
+{
+  // Do not handle transformation boundary conditions on lower universes
+  if (n_coord() != 1) {
+    mark_as_lost("Cannot transform particle " + std::to_string(id()) +
+                 " off surface in a lower universe.");
+    return;
+  }
+
+  // Score surface currents since transformation causes the direction of the
+  // particle to change. For surface filters, we need to score the tallies
+  // twice, once before the particle's surface attribute has changed and
+  // once after. For mesh surface filters, we need to artificially move
+  // the particle slightly back in case the surface crossing is coincident
+  // with a mesh boundary
+  if (!model::active_surface_tallies.empty()) {
+    score_surface_tally(*this, model::active_surface_tallies);
+  }
+
+  if (!model::active_meshsurf_tallies.empty()) {
+    Position r {this->r()};
+    this->r() -= TINY_BIT * u();
+    score_surface_tally(*this, model::active_meshsurf_tallies);
+    this->r() = r;
+  }
+
+  // Adjust the particle's location and direction.
+  r() = new_r;
+  u() = new_u;
+
+  // Clear the surface assignment after transformation so as not to confuse the
+  // cell finding routine
+  surface() = SURFACE_NONE;
+
+  // Figure out what cell particle is in now
+  // If a transformation surface is coincident with a lattice or universe
+  // boundary, it is necessary to redetermine the particle's coordinates in
+  // the lower universes.
+  // (unless we're using a dagmc model, which has exactly one universe)
+  n_coord() = 1;
+  if (surf.geom_type() != GeometryType::DAG && !exhaustive_find_cell(*this)) {
+    mark_as_lost("Couldn't find particle after transforming from surface " +
+                 std::to_string(surf.id_) + ".");
+    return;
+  }
+
+  // Set previous coordinate going slightly past surface crossing
+  r_last_current() = r() + TINY_BIT * u();
+
+  // Diagnostic message
+  if (settings::verbosity >= 10 || trace()) {
+    write_message(1, "    Transformed from surface {}", surf.id_);
+  }
+}
+
 void Particle::mark_as_lost(const char* message)
 {
   // Print warning and write lost particle file

src/surface.cpp

@@ -88,6 +88,39 @@ Surface::Surface(pugi::xml_node surf_node)
       bc_ = make_unique<WhiteBC>();
     } else if (surf_bc == "periodic") {
       // Periodic BCs are handled separately
+    } else if (surf_bc == "transformation") {
+      int vector_size_exp = 12;
+      vector<double> dir_trans(vector_size_exp);
+      vector<double> pos_trans(vector_size_exp);
+
+      if (check_for_node(surf_node, "direction_transformation")) {
+        dir_trans =
+          get_node_array<double>(surf_node, "direction_transformation");
+        if (dir_trans.size() != vector_size_exp) {
+          fatal_error(fmt::format(
+            "Transformation on surface {} expects direction matrix size {} "
+            "but was given size {}",
+            id_, vector_size_exp, dir_trans.size()));
+        }
+      } else {
+        dir_trans = {
+          1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0};
+      }
+
+      if (check_for_node(surf_node, "position_transformation")) {
+        pos_trans =
+          get_node_array<double>(surf_node, "position_transformation");
+        if (pos_trans.size() != vector_size_exp) {
+          fatal_error(
+            fmt::format("Transformation on surface {} expects position matrix "
+                        "size {} but was given size {}",
+              id_, vector_size_exp, pos_trans.size()));
+        }
+      } else {
+        pos_trans = {
+          1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0};
+      }
+      bc_ = make_unique<TransformationBC>(dir_trans, pos_trans);
     } else {
       fatal_error(fmt::format("Unknown boundary condition \"{}\" specified "
                               "on surface {}",

tests/regression_tests/transformation_direction/inputs_true.dat

@@ -0,0 +1,29 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" depletable="true">
+      <density value="4.5" units="g/cc"/>
+      <nuclide name="U235" ao="1.0"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="1 -2 3 -4 5 -6" universe="1"/>
+    <surface id="1" type="x-plane" boundary="transformation" coeffs="-5.0" direction_transformation="-1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0" position_transformation="1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0"/>
+    <surface id="2" type="x-plane" boundary="transformation" coeffs="5.0" direction_transformation="-1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0" position_transformation="1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0"/>
+    <surface id="3" type="y-plane" boundary="transformation" coeffs="-5.0" direction_transformation="1.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 1.0 0.0" position_transformation="1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0"/>
+    <surface id="4" type="y-plane" boundary="transformation" coeffs="5.0" direction_transformation="1.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 1.0 0.0" position_transformation="1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0"/>
+    <surface id="5" type="z-plane" boundary="transformation" coeffs="-5.0" direction_transformation="1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 -1.0 0.0" position_transformation="1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0"/>
+    <surface id="6" type="z-plane" boundary="transformation" coeffs="5.0" direction_transformation="1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 -1.0 0.0" position_transformation="1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <source type="independent" strength="1.0" particle="neutron">
+      <space type="box">
+        <parameters>-4.0 -4.0 -4.0 4.0 4.0 4.0</parameters>
+      </space>
+    </source>
+  </settings>
+  </model>

tests/regression_tests/transformation_direction/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+2.279066E+00 4.793565E-03

tests/regression_tests/transformation_direction/test.py

@@ -0,0 +1,92 @@
+import openmc
+import pytest
+import numpy as np
+
+from tests.testing_harness import PyAPITestHarness
+
+
+@pytest.fixture
+def box_model():
+    model = openmc.model.Model()
+
+    fuel = openmc.Material()
+    fuel.add_nuclide('U235', 1.0)
+    fuel.set_density('g/cc', 4.5)
+
+    normal_x = np.array([1.0, 0.0, 0.0])
+    dir_transform_xplane = np.identity(3) - 2 * np.outer(
+        normal_x, normal_x
+    )
+    dir_transform_xplane = np.append(dir_transform_xplane, np.zeros((3,1)), axis=1).flatten()
+
+    normal_y = np.array([0.0, 1.0, 0.0])
+    dir_transform_yplane = np.identity(3) - 2 * np.outer(
+        normal_y, normal_y
+    )
+    dir_transform_yplane = np.append(dir_transform_yplane, np.zeros((3,1)), axis=1).flatten()
+
+    normal_z = np.array([0.0, 0.0, 1.0])
+    dir_transform_zplane = np.identity(3) - 2 * np.outer(
+        normal_z, normal_z
+    )
+    dir_transform_zplane = np.append(dir_transform_zplane, np.zeros((3,1)), axis=1).flatten()
+
+    neg_xplane = openmc.XPlane(
+        x0=-5.0,
+        boundary_type="transformation",
+        transformation={"direction": dir_transform_xplane}
+    )
+    pos_xplane = openmc.XPlane(
+        x0=5.0,
+        boundary_type="transformation",
+        transformation={"direction": dir_transform_xplane}
+    )
+
+    neg_yplane = openmc.YPlane(
+        y0=-5.0,
+        boundary_type="transformation",
+        transformation={"direction": dir_transform_yplane}
+    )
+    pos_yplane = openmc.YPlane(
+        y0=5.0,
+        boundary_type="transformation",
+        transformation={"direction": dir_transform_yplane}
+    )
+
+    neg_zplane = openmc.ZPlane(
+        z0=-5.0,
+        boundary_type="transformation",
+        transformation={"direction": dir_transform_zplane}
+    )
+    pos_zplane = openmc.ZPlane(
+        z0=5.0,
+        boundary_type="transformation",
+        transformation={"direction": dir_transform_zplane}
+    )
+
+    region = (
+        +neg_xplane
+        & -pos_xplane
+        & +neg_yplane
+        & -pos_yplane
+        & +neg_zplane
+        & -pos_zplane
+    )
+    cell = openmc.Cell(fill=fuel, region=region)
+
+    model.geometry = openmc.Geometry([cell])
+
+    model.settings.particles = 1000
+    model.settings.batches = 10
+    model.settings.inactive = 5
+
+    source = openmc.IndependentSource()
+    source.space = openmc.stats.Box([-4.0, -4.0, -4.0], [4.0, 4.0, 4.0])
+    model.settings.source = source
+
+    return model
+
+
+def test_transformation_direction(box_model):
+    harness = PyAPITestHarness('statepoint.10.h5', box_model)
+    harness.main()

tests/regression_tests/transformation_mixed/inputs_true.dat

@@ -0,0 +1,37 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1">
+      <density value="1.0" units="g/cc"/>
+      <nuclide name="H1" ao="2.0"/>
+      <nuclide name="O16" ao="1.0"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+    <material id="2" depletable="true">
+      <density value="4.5" units="g/cc"/>
+      <nuclide name="U235" ao="1.0"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="1 -2 3 -4 5 -6 7" universe="0"/>
+    <cell id="2" material="2" region="3 5 -6 -7" universe="0"/>
+    <surface id="1" type="x-plane" boundary="transformation" coeffs="0.0" direction_transformation="0.0 1.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0" position_transformation="0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0"/>
+    <surface id="2" type="x-plane" boundary="transformation" coeffs="5.0" direction_transformation="-1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0" position_transformation="1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0"/>
+    <surface id="3" type="y-plane" boundary="transformation" coeffs="0.0" direction_transformation="0.0 -1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0" position_transformation="0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0"/>
+    <surface id="4" type="y-plane" boundary="transformation" coeffs="5.0" direction_transformation="1.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 1.0 0.0" position_transformation="1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0"/>
+    <surface id="5" type="z-plane" boundary="transformation" coeffs="-5.0" direction_transformation="1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0" position_transformation="1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 10.0"/>
+    <surface id="6" type="plane" boundary="transformation" coeffs="0 0 1 5.0" direction_transformation="1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0" position_transformation="1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 -10.0"/>
+    <surface id="7" type="z-cylinder" coeffs="2.5 0.0 2.0"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>4</batches>
+    <inactive>0</inactive>
+    <source type="independent" strength="1.0" particle="neutron">
+      <space type="box">
+        <parameters>0 0 0 5 5 0</parameters>
+      </space>
+    </source>
+  </settings>
+  </model>

tests/regression_tests/transformation_mixed/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+1.624889E+00 1.108153E-02