diff --git a/doc/source/userguide/post_processors.rst b/doc/source/userguide/post_processors.rst
index ecefe1491..70d42f115 100644
--- a/doc/source/userguide/post_processors.rst
+++ b/doc/source/userguide/post_processors.rst
@@ -535,6 +535,131 @@ With energy restriction:
# result[0] is a vector with one element (the single group)
print(result[0][0])
+Surface Postprocessor
+=====================
+
+The :py:class:`pyopensn.post.SurfacePostprocessor` computes surface flux/current
+integrals, maxima, minima, or area-weighted averages over boundary surfaces and
+energy groups. It is specifically designed for discrete ordinates problems
+and provides access to boundary partial currents and net currents.
+
+Basic Usage
+-----------
+
+Create and execute a surface postprocessor:
+
+.. code-block:: python
+
+ from pyopensn.post import SurfacePostprocessor
+
+ pps = SurfacePostprocessor(
+ problem=phys,
+ value_type="integral",
+ current_type="net",
+ boundaries=["xmin"],
+ )
+ pps.Execute()
+ result = pps.GetValue()
+
+Operation Types
+---------------
+
+The ``value_type`` parameter determines the reduction operation:
+
+* ``"integral"`` — area-weighted integral of the selected current type
+* ``"avg"`` — area-weighted average of the selected current type
+* ``"max"`` — maximum value of the selected current type on the surface
+* ``"min"`` — minimum value of the selected current type on the surface
+
+The ``current_type`` parameter selects which component of the surface flux is
+evaluated:
+
+* ``"incoming"`` — partial current entering the domain
+* ``"outgoing"`` — partial current leaving the domain
+* ``"net"`` — net current (outgoing minus incoming)
+
+Spatial Restriction
+-------------------
+
+A surface postprocessor must be restricted to one or more boundaries:
+
+.. code-block:: python
+
+ pps = SurfacePostprocessor(
+ problem=phys,
+ value_type="integral",
+ current_type="net",
+ boundaries=["xmin", "xmax"],
+ )
+
+Optionally, restrict computation to one or more logical volumes to evaluate
+only a portion of the specified boundaries:
+
+.. code-block:: python
+
+ from pyopensn.logvol import RPPLogicalVolume
+
+ lv = RPPLogicalVolume(ymin=0.0, ymax=0.5, infx=True, infz=True)
+ pps = SurfacePostprocessor(
+ problem=phys,
+ value_type="integral",
+ current_type="net",
+ boundaries=["xmin"],
+ logical_volumes=[lv],
+ )
+
+Each logical volume produces one row of results. If no logical volumes are
+specified, the entire surface of the named boundaries is used as a single region.
+
+Energy Restriction
+-------------------
+
+Similar to the volume postprocessor, you can restrict the computation to a
+single group or a single groupset:
+
+.. code-block:: python
+
+ pps = SurfacePostprocessor(
+ problem=phys,
+ value_type="integral",
+ current_type="net",
+ boundaries=["xmin"],
+ group=6,
+ )
+
+Or:
+
+.. code-block:: python
+
+ pps = SurfacePostprocessor(
+ problem=phys,
+ value_type="integral",
+ current_type="net",
+ boundaries=["xmin"],
+ groupset=1,
+ )
+
+Results
+-------
+
+After calling ``Execute()``, retrieve results with ``GetValue()``. The return
+value is a 2D array indexed as ``result[region][group]``:
+
+.. code-block:: python
+
+ pps = SurfacePostprocessor(
+ problem=phys,
+ value_type="integral",
+ current_type="net",
+ boundaries=["xmin"],
+ )
+ pps.Execute()
+ result = pps.GetValue()
+
+ # result[0] is a vector of values for the combined xmin boundary, one per group
+ for group_index, value in enumerate(result[0]):
+ print(f"Group {group_index}: {value}")
+
Other Useful Post-Processing Paths
==================================
diff --git a/modules/linear_boltzmann_solvers/lbs_problem/postprocessors/surface_postprocessor.cc b/modules/linear_boltzmann_solvers/lbs_problem/postprocessors/surface_postprocessor.cc
new file mode 100644
index 000000000..10e4a7782
--- /dev/null
+++ b/modules/linear_boltzmann_solvers/lbs_problem/postprocessors/surface_postprocessor.cc
@@ -0,0 +1,455 @@
+// SPDX-FileCopyrightText: 2026 The OpenSn Authors
+// SPDX-License-Identifier: MIT
+
+#include "modules/linear_boltzmann_solvers/lbs_problem/postprocessors/surface_postprocessor.h"
+#include "framework/object_factory.h"
+#include "framework/math/spatial_discretization/finite_element/finite_element_data.h"
+#include "framework/runtime.h"
+#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.h"
+#include
+#include
+#include
+#include
+
+namespace opensn
+{
+
+namespace
+{
+
+auto
+ComputeInflow(const LBSGroupset& groupset,
+ const Cell& cell,
+ const CellMapping& cell_mapping,
+ SweepBoundary& bndry,
+ const std::vector>& IntS_shapeI,
+ const CellFace& face,
+ unsigned int face_idx,
+ unsigned int g)
+{
+ double inflow = 0;
+ const auto& omegas = groupset.quadrature->GetOmegas();
+ for (int n = 0; n < omegas.size(); ++n)
+ {
+ const double wt = groupset.quadrature->GetWeight(n);
+ const double mu = omegas[n].Dot(face.normal);
+ if (mu >= 0.)
+ continue;
+
+ for (int fi = 0; fi < face.vertex_ids.size(); ++fi)
+ {
+ const int i = cell_mapping.MapFaceNode(face_idx, fi);
+ const auto& IntFi_shapeI = IntS_shapeI[face_idx](i);
+
+ const double* psi_ptr = bndry.PsiIncoming(cell.local_id, face_idx, fi, n, groupset.id, g);
+ const double psi = (psi_ptr != nullptr) ? *psi_ptr : 0.0;
+ inflow += mu * wt * psi * IntFi_shapeI;
+ }
+ }
+ return inflow;
+};
+} // namespace
+
+OpenSnRegisterObjectInNamespace(lbs, SurfacePostprocessor);
+
+InputParameters
+SurfacePostprocessor::GetInputParameters()
+{
+ InputParameters params;
+ params.AddRequiredParameter>("problem",
+ "A handle to an existing LBS problem.");
+ params.AddRequiredParameterArray(
+ "boundaries",
+ "Boundary restriction for the postprocessor. Empty/unspecified means no block restriction.");
+ params.AddOptionalParameterArray("logical_volumes",
+ std::vector>{},
+ "Logical volume to restrict the computation to.");
+ params.AddRequiredParameter(
+ "value_type", "Type of value to compute: 'integral', 'max', 'min', or 'avg'");
+ params.AddRequiredParameter(
+ "current_type", "Type of value to compute: 'incoming', 'outgoing', 'net'");
+ params.AddOptionalParameter(
+ "group", 0, "Single group to compute (mutually exclusive with groupset).");
+ params.AddOptionalParameter(
+ "groupset", 0, "Single groupset to compute (mutually exclusive with group).");
+
+ params.ConstrainParameterRange("value_type",
+ AllowableRangeList::New({"integral", "max", "min", "avg"}));
+ params.ConstrainParameterRange("current_type",
+ AllowableRangeList::New({"incoming", "outgoing", "net"}));
+
+ return params;
+}
+
+std::shared_ptr
+SurfacePostprocessor::Create(const ParameterBlock& params)
+{
+ auto& factory = opensn::ObjectFactory::GetInstance();
+ return factory.Create("lbs::SurfacePostprocessor", params);
+}
+
+SurfacePostprocessor::SurfacePostprocessor(const InputParameters& params)
+ : do_problem_(params.GetSharedPtrParam("problem")),
+ boundary_names_(params.GetParamVectorValue("boundaries")),
+ logical_volumes_(params.GetParamVectorValue>("logical_volumes")),
+ selected_group_(params.IsParameterValid("group")
+ ? std::make_optional(params.GetParamValue("group"))
+ : std::nullopt),
+ selected_groupset_(params.IsParameterValid("groupset")
+ ? std::make_optional(params.GetParamValue("groupset"))
+ : std::nullopt)
+{
+ const auto& grid = do_problem_->GetGrid();
+ const auto& bnd_id_map = grid->GetBoundaryNameMap();
+ for (auto& name : boundary_names_)
+ {
+ const auto& id = bnd_id_map.at(name);
+ boundary_ids_.push_back(id);
+ }
+ // const auto & bnd_ids = do_problem_->GetSweepBoundaries();
+
+ if (selected_group_.has_value() && selected_groupset_.has_value())
+ throw std::invalid_argument("'group' and 'groupset' cannot be specified together");
+
+ if (selected_group_.has_value() && selected_group_.value() >= do_problem_->GetNumGroups())
+ throw std::invalid_argument("'group' must be less than " +
+ std::to_string(do_problem_->GetNumGroups()));
+
+ if (selected_groupset_.has_value() &&
+ selected_groupset_.value() >= do_problem_->GetGroupsets().size())
+ throw std::invalid_argument("'groupset' must be less than " +
+ std::to_string(do_problem_->GetGroupsets().size()));
+
+ const auto value_type_str = params.GetParamValue("value_type");
+ if (value_type_str == "max")
+ value_type_ = ValueType::MAX;
+ else if (value_type_str == "min")
+ value_type_ = ValueType::MIN;
+ else if (value_type_str == "integral")
+ value_type_ = ValueType::INTEGRAL;
+ else if (value_type_str == "avg")
+ value_type_ = ValueType::AVERAGE;
+
+ const auto current_type_str = params.GetParamValue("current_type");
+ if (current_type_str == "incoming")
+ current_type_ = CurrentType::INCOMING;
+ else if (current_type_str == "outgoing")
+ current_type_ = CurrentType::OUTGOING;
+ else if (current_type_str == "net")
+ current_type_ = CurrentType::NET;
+
+ CreateSpatialRestriction();
+ CreateEnergyRestriction();
+
+ auto n_lvs = std::max(static_cast(1), logical_volumes_.size());
+ values_.resize(n_lvs);
+ for (auto& vals : values_)
+ vals.resize(groups_.size());
+}
+
+void
+SurfacePostprocessor::CreateSpatialRestriction()
+{
+ const auto& grid = do_problem_->GetGrid();
+ const auto& groupsets = do_problem_->GetGroupsets();
+
+ std::vector all_side_faces;
+ for (const auto& cell : grid->local_cells)
+ {
+ for (int f = 0; f < cell.faces.size(); ++f)
+ {
+ const auto& face = cell.faces[f];
+ // boundary face?
+ if (not face.has_neighbor)
+ {
+ if (std::find(boundary_ids_.begin(), boundary_ids_.end(), face.neighbor_id) !=
+ boundary_ids_.end())
+ all_side_faces.push_back({cell.local_id, f});
+ }
+ }
+ }
+
+ if (logical_volumes_.empty())
+ {
+ side_faces_.resize(1);
+ side_faces_[0] = all_side_faces;
+ }
+ else
+ {
+ side_faces_.resize(logical_volumes_.size());
+ for (unsigned int i = 0; i < logical_volumes_.size(); ++i)
+ side_faces_[i] = GetLogicalVolumeSides(all_side_faces, logical_volumes_[i]);
+ }
+}
+
+std::vector
+SurfacePostprocessor::GetLogicalVolumeSides(const std::vector& all_side_faces,
+ std::shared_ptr log_vol)
+{
+ const auto& grid = do_problem_->GetGrid();
+ std::vector side_faces;
+ for (const auto& [cell, face_idx] : all_side_faces)
+ {
+ const auto& face = grid->local_cells[cell].faces[face_idx];
+ if (log_vol->Inside(face.centroid))
+ side_faces.emplace_back(cell, face_idx);
+ }
+ return side_faces;
+}
+
+void
+SurfacePostprocessor::CreateEnergyRestriction()
+{
+ if (selected_group_.has_value())
+ {
+ const auto g = selected_group_.value();
+ groups_.push_back(g);
+ for (unsigned int gs_id = 0; gs_id < do_problem_->GetNumGroupsets(); gs_id++)
+ {
+ const auto& gs = do_problem_->GetGroupset(gs_id);
+ if ((gs.first_group <= g) and (g <= gs.last_group))
+ {
+ groupset_ids_.push_back(gs_id);
+ break;
+ }
+ }
+ }
+ else if (selected_groupset_.has_value())
+ {
+ const auto groupset_id = selected_groupset_.value();
+ const auto& groupset = do_problem_->GetGroupsets()[groupset_id];
+ for (unsigned int g = groupset.first_group; g <= groupset.last_group; ++g)
+ {
+ groups_.push_back(g);
+ groupset_ids_.push_back(groupset_id);
+ }
+ }
+ else
+ {
+ const auto& groupsets = do_problem_->GetGroupsets();
+ for (unsigned int gsi = 0; gsi < groupsets.size(); ++gsi)
+ {
+ const auto& gs = groupsets[gsi];
+ for (unsigned int g = gs.first_group; g <= gs.last_group; g++)
+ {
+ groups_.push_back(g);
+ groupset_ids_.push_back(gsi);
+ }
+ }
+ }
+}
+
+void
+SurfacePostprocessor::Execute()
+{
+ for (unsigned int i = 0; i < side_faces_.size(); ++i)
+ {
+ switch (value_type_)
+ {
+ case ValueType::INTEGRAL:
+ values_[i] = ComputeIntegral(side_faces_[i]);
+ break;
+
+ case ValueType::MAX:
+ values_[i] = ComputeMax(side_faces_[i]);
+ break;
+
+ case ValueType::MIN:
+ values_[i] = ComputeMin(side_faces_[i]);
+ break;
+
+ case ValueType::AVERAGE:
+ values_[i] = ComputeAvg(side_faces_[i]);
+ break;
+ }
+ }
+}
+
+std::vector
+SurfacePostprocessor::ComputeIntegral(const std::vector& side_faces)
+{
+ const auto& sdm = do_problem_->GetSpatialDiscretization();
+ const auto& grid = sdm.GetGrid();
+ const auto& sweep_boundaries = do_problem_->GetSweepBoundaries();
+ const auto& cell_outflow_views = do_problem_->GetCellOutflowViews();
+ const auto& unit_cell_matrices = do_problem_->GetUnitCellMatrices();
+ const auto& groupsets = do_problem_->GetGroupsets();
+ const auto& uk_man = do_problem_->GetUnknownManager();
+ const auto phi = do_problem_->GetPhiNewLocal();
+ auto coord = sdm.GetSpatialWeightingFunction();
+
+ std::vector local_integral(groups_.size(), 0.0);
+ for (const auto [cell_id, face_idx] : side_faces)
+ {
+ const auto& cell = grid->local_cells[cell_id];
+ const auto& cell_mapping = sdm.GetCellMapping(cell);
+ const auto& outflow_view = cell_outflow_views[cell.local_id];
+ const auto& fe_intgrl_values = unit_cell_matrices[cell.local_id];
+ const auto& IntS_shapeI = fe_intgrl_values.intS_shapeI;
+
+ const auto& face = cell.faces[face_idx];
+ const auto& bndry = sweep_boundaries.at(face.neighbor_id);
+ for (std::size_t k = 0; k < groups_.size(); ++k)
+ {
+ const auto g = groups_[k];
+ const auto& groupset = groupsets[groupset_ids_[k]];
+
+ if (current_type_ == CurrentType::INCOMING)
+ local_integral[k] +=
+ ComputeInflow(groupset, cell, cell_mapping, *bndry, IntS_shapeI, face, face_idx, g);
+ else if (current_type_ == CurrentType::OUTGOING)
+ local_integral[k] += outflow_view.Get(face_idx, g);
+ else if (current_type_ == CurrentType::NET)
+ local_integral[k] +=
+ outflow_view.Get(face_idx, g) +
+ ComputeInflow(groupset, cell, cell_mapping, *bndry, IntS_shapeI, face, face_idx, g);
+ }
+ }
+
+ std::vector global_integral(groups_.size(), 0.0);
+ for (std::size_t i = 0; i < local_integral.size(); ++i)
+ mpi_comm.all_reduce(local_integral[i], global_integral[i], mpi::op::sum());
+
+ return global_integral;
+}
+
+std::vector
+SurfacePostprocessor::ComputeMin(const std::vector& side_faces)
+{
+ const auto& sdm = do_problem_->GetSpatialDiscretization();
+ const auto& grid = sdm.GetGrid();
+ const auto& sweep_boundaries = do_problem_->GetSweepBoundaries();
+ const auto& cell_outflow_views = do_problem_->GetCellOutflowViews();
+ const auto& unit_cell_matrices = do_problem_->GetUnitCellMatrices();
+ const auto& groupsets = do_problem_->GetGroupsets();
+ const auto& uk_man = do_problem_->GetUnknownManager();
+ const auto phi = do_problem_->GetPhiNewLocal();
+ auto coord = sdm.GetSpatialWeightingFunction();
+
+ std::vector local_values(groups_.size(), std::numeric_limits::infinity());
+ for (const auto [cell_id, face_idx] : side_faces)
+ {
+ const auto& cell = grid->local_cells[cell_id];
+ const auto& cell_mapping = sdm.GetCellMapping(cell);
+ const auto& outflow_view = cell_outflow_views[cell.local_id];
+ const auto& fe_intgrl_values = unit_cell_matrices[cell.local_id];
+ const auto& IntS_shapeI = fe_intgrl_values.intS_shapeI;
+
+ const auto& face = cell.faces[face_idx];
+ const auto& bndry = sweep_boundaries.at(face.neighbor_id);
+ for (std::size_t k = 0; k < groups_.size(); ++k)
+ {
+ const auto g = groups_[k];
+ const auto& groupset = groupsets[groupset_ids_[k]];
+
+ if (current_type_ == CurrentType::INCOMING)
+ local_values[k] = std::min(
+ local_values[k],
+ ComputeInflow(groupset, cell, cell_mapping, *bndry, IntS_shapeI, face, face_idx, g));
+ else if (current_type_ == CurrentType::OUTGOING)
+ local_values[k] = std::min(local_values[k], outflow_view.Get(face_idx, g));
+ else if (current_type_ == CurrentType::NET)
+ local_values[k] = std::min(
+ local_values[k],
+ outflow_view.Get(face_idx, g) +
+ ComputeInflow(groupset, cell, cell_mapping, *bndry, IntS_shapeI, face, face_idx, g));
+ }
+ }
+
+ std::vector global_values(groups_.size(), 0.0);
+ for (std::size_t i = 0; i < local_values.size(); ++i)
+ mpi_comm.all_reduce(local_values[i], global_values[i], mpi::op::min());
+
+ return global_values;
+}
+
+std::vector
+SurfacePostprocessor::ComputeMax(const std::vector& side_faces)
+{
+ const auto& sdm = do_problem_->GetSpatialDiscretization();
+ const auto& grid = sdm.GetGrid();
+ const auto& sweep_boundaries = do_problem_->GetSweepBoundaries();
+ const auto& cell_outflow_views = do_problem_->GetCellOutflowViews();
+ const auto& unit_cell_matrices = do_problem_->GetUnitCellMatrices();
+ const auto& groupsets = do_problem_->GetGroupsets();
+ const auto& uk_man = do_problem_->GetUnknownManager();
+ const auto phi = do_problem_->GetPhiNewLocal();
+ auto coord = sdm.GetSpatialWeightingFunction();
+
+ std::vector local_values(groups_.size(), -std::numeric_limits::infinity());
+ for (const auto [cell_id, face_idx] : side_faces)
+ {
+ const auto& cell = grid->local_cells[cell_id];
+ const auto& cell_mapping = sdm.GetCellMapping(cell);
+ const auto& outflow_view = cell_outflow_views[cell.local_id];
+ const auto& fe_intgrl_values = unit_cell_matrices[cell.local_id];
+ const auto& IntS_shapeI = fe_intgrl_values.intS_shapeI;
+
+ const auto& face = cell.faces[face_idx];
+ const auto& bndry = sweep_boundaries.at(face.neighbor_id);
+ for (std::size_t k = 0; k < groups_.size(); ++k)
+ {
+ const auto g = groups_[k];
+ const auto& groupset = groupsets[groupset_ids_[k]];
+
+ if (current_type_ == CurrentType::INCOMING)
+ local_values[k] = std::max(
+ local_values[k],
+ ComputeInflow(groupset, cell, cell_mapping, *bndry, IntS_shapeI, face, face_idx, g));
+ else if (current_type_ == CurrentType::OUTGOING)
+ local_values[k] = std::max(local_values[k], outflow_view.Get(face_idx, g));
+ else if (current_type_ == CurrentType::NET)
+ local_values[k] = std::max(
+ local_values[k],
+ outflow_view.Get(face_idx, g) +
+ ComputeInflow(groupset, cell, cell_mapping, *bndry, IntS_shapeI, face, face_idx, g));
+ }
+ }
+
+ std::vector global_values(groups_.size(), 0.0);
+ for (std::size_t i = 0; i < local_values.size(); ++i)
+ mpi_comm.all_reduce(local_values[i], global_values[i], mpi::op::max());
+
+ return global_values;
+}
+
+std::vector
+SurfacePostprocessor::ComputeAvg(const std::vector& side_faces)
+{
+ const auto& sdm = do_problem_->GetSpatialDiscretization();
+ const auto& grid = sdm.GetGrid();
+ auto coord = sdm.GetSpatialWeightingFunction();
+
+ double local_weighted_area = 0.0;
+ for (const auto [cell_id, face_idx] : side_faces)
+ {
+ const auto& cell = grid->local_cells[cell_id];
+ const auto& cell_mapping = sdm.GetCellMapping(cell);
+ const auto fe_face_data = cell_mapping.MakeSurfaceFiniteElementData(face_idx);
+
+ for (const std::size_t qp : fe_face_data.GetQuadraturePointIndices())
+ local_weighted_area += coord(fe_face_data.QPointXYZ(qp)) * fe_face_data.JxW(qp);
+ }
+ double global_weighted_area = 0.0;
+ mpi_comm.all_reduce(local_weighted_area, global_weighted_area, mpi::op::sum());
+
+ auto global_weighted_integral = ComputeIntegral(side_faces);
+
+ std::vector values(groups_.size());
+ for (std::size_t i = 0; i < global_weighted_integral.size(); ++i)
+ {
+ if (global_weighted_area > 0.0)
+ values[i] = global_weighted_integral[i] / global_weighted_area;
+ else
+ values[i] = 0.0;
+ }
+ return values;
+}
+
+std::vector>
+SurfacePostprocessor::GetValue() const
+{
+ return values_;
+}
+
+} // namespace opensn
diff --git a/modules/linear_boltzmann_solvers/lbs_problem/postprocessors/surface_postprocessor.h b/modules/linear_boltzmann_solvers/lbs_problem/postprocessors/surface_postprocessor.h
new file mode 100644
index 000000000..41031d925
--- /dev/null
+++ b/modules/linear_boltzmann_solvers/lbs_problem/postprocessors/surface_postprocessor.h
@@ -0,0 +1,86 @@
+// SPDX-FileCopyrightText: 2026 The OpenSn Authors
+// SPDX-License-Identifier: MIT
+
+#pragma once
+
+#include "framework/parameters/input_parameters.h"
+#include "framework/mesh/logical_volume/logical_volume.h"
+#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.h"
+#include
+
+namespace opensn
+{
+
+class SurfacePostprocessor
+{
+public:
+ enum class ValueType
+ {
+ INTEGRAL,
+ MAX,
+ MIN,
+ AVERAGE
+ };
+
+ enum class CurrentType
+ {
+ NET,
+ INCOMING,
+ OUTGOING
+ };
+
+ /// Input parameters based construction.
+ explicit SurfacePostprocessor(const InputParameters& params);
+
+ void Execute();
+
+ std::vector> GetValue() const;
+
+private:
+ struct SideFace
+ {
+ uint64_t cell_local_index;
+ int local_face;
+ };
+
+ void CreateSpatialRestriction();
+ void CreateEnergyRestriction();
+ std::vector
+ GetLogicalVolumeSides(const std::vector& all_side_faces,
+ std::shared_ptr log_vol);
+ std::vector ComputeIntegral(const std::vector& side_faces);
+ std::vector ComputeMin(const std::vector& side_faces);
+ std::vector ComputeMax(const std::vector& side_faces);
+ std::vector ComputeAvg(const std::vector& side_faces);
+
+ std::shared_ptr do_problem_;
+ ///
+ std::vector boundary_names_;
+ /// Boundary IDs this postprocessor is restricted to
+ std::vector boundary_ids_;
+ /// Logical volume associated with this PPS (can be null)
+ std::vector> logical_volumes_;
+ ///
+ std::vector> side_faces_;
+ /// Groups
+ std::vector groups_;
+ /// Groupset IDs
+ std::vector groupset_ids_;
+ /// Type of value to compute
+ ValueType value_type_;
+ /// Type of current
+ CurrentType current_type_;
+ /// Computed postprocessed values
+ std::vector> values_;
+ /// Selected group (-1 means not selected)
+ std::optional selected_group_;
+ /// Selected groupset (-1 means not selected)
+ std::optional selected_groupset_;
+
+public:
+ /// Returns the input parameters for this object.
+ static InputParameters GetInputParameters();
+ static std::shared_ptr Create(const ParameterBlock& params);
+};
+
+} // namespace opensn
diff --git a/python/lib/post.cc b/python/lib/post.cc
index 84082549f..cdec0b09d 100644
--- a/python/lib/post.cc
+++ b/python/lib/post.cc
@@ -3,6 +3,7 @@
#include "python/lib/py_wrappers.h"
#include "modules/linear_boltzmann_solvers/lbs_problem/postprocessors/volume_postprocessor.h"
+#include "modules/linear_boltzmann_solvers/lbs_problem/postprocessors/surface_postprocessor.h"
#include
#include
@@ -73,6 +74,72 @@ WrapPostprocessors(py::module& post)
Columns correspond to the energy restrictions (i.e. groups, or groups within a groupset if specified)
)"
);
+
+ // Surface post-processor value type enum
+ py::enum_(post, "SurfacePostprocessorValueType")
+ .value("INTEGRAL", SurfacePostprocessor::ValueType::INTEGRAL)
+ .value("MAX", SurfacePostprocessor::ValueType::MAX)
+ .value("MIN", SurfacePostprocessor::ValueType::MIN)
+ .value("AVERAGE", SurfacePostprocessor::ValueType::AVERAGE);
+
+ py::enum_(post, "SurfacePostprocessorCurrentType")
+ .value("INCOMING", SurfacePostprocessor::CurrentType::INCOMING)
+ .value("OUTGOING", SurfacePostprocessor::CurrentType::OUTGOING)
+ .value("NET", SurfacePostprocessor::CurrentType::NET);
+
+ // Surface post-processor
+ auto sp = py::class_>(
+ post,
+ "SurfacePostprocessor",
+ R"(
+ Surface post-processor.
+
+ Wrapper of :cpp:class:`opensn::SurfacePostprocessor`.
+ )"
+ );
+ sp.def(
+ py::init(
+ [](py::kwargs& params)
+ {
+ return SurfacePostprocessor::Create(kwargs_to_param_block(params));
+ }
+ ),
+ R"(
+ Construct a surface post-processor object.
+
+ Parameters
+ ----------
+ problem : DiscreteOrdinatesProblem
+ A handle to an existing discrete ordinates problem.
+ value_type : str, required
+ Type of value to compute: 'integral', 'max', 'min', or 'avg'.
+ current_type : str, required
+ Type of value to compute: 'incoming', 'outgoing', or 'net'.
+ )"
+ );
+ sp.def(
+ "Execute",
+ [](SurfacePostprocessor& self){
+ self.Execute();
+ },
+ R"(
+ Execute the postprocessor
+ )"
+ );
+ sp.def(
+ "GetValue",
+ [](SurfacePostprocessor& self)
+ {
+ return self.GetValue();
+ },
+ R"(
+ Returns the value of the postprocessor.
+
+ Rows correspond to the spatial restriction (i.e. logical volumes, if specified)
+ Columns correspond to the energy restrictions (i.e. groups, or groups within a groupset if specified)
+ )"
+ );
+
// clang-format on
}
diff --git a/test/python/modules/linear_boltzmann_solvers/pps/tests.json b/test/python/modules/linear_boltzmann_solvers/pps/tests.json
index 774af92f9..7274dc227 100644
--- a/test/python/modules/linear_boltzmann_solvers/pps/tests.json
+++ b/test/python/modules/linear_boltzmann_solvers/pps/tests.json
@@ -22,5 +22,17 @@
"error_code": 0
}
]
+ },
+ {
+ "file": "transport_2d_1_poly.py",
+ "comment": "Computing surface post-processors on 2D domain",
+ "num_procs": 4,
+ "weight_class": "intermediate",
+ "checks": [
+ {
+ "type": "ErrorCode",
+ "error_code": 0
+ }
+ ]
}
]
diff --git a/test/python/modules/linear_boltzmann_solvers/pps/transport_2d_1_poly.py b/test/python/modules/linear_boltzmann_solvers/pps/transport_2d_1_poly.py
new file mode 100644
index 000000000..92879d08b
--- /dev/null
+++ b/test/python/modules/linear_boltzmann_solvers/pps/transport_2d_1_poly.py
@@ -0,0 +1,213 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+"""
+2D PWLD transport test with vacuum and Incident-isotropic bundary conditions
+Test: Max-value=0.50758 and 2.52527e-04
+"""
+
+import math
+import os
+import sys
+
+if "opensn_console" not in globals():
+ from mpi4py import MPI
+
+ size = MPI.COMM_WORLD.size
+ rank = MPI.COMM_WORLD.rank
+ sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), "../../../../../")))
+ from pyopensn.aquad import GLCProductQuadrature2DXY
+ from pyopensn.logvol import RPPLogicalVolume
+ from pyopensn.mesh import FromFileMeshGenerator, KBAGraphPartitioner
+ from pyopensn.post import SurfacePostprocessor
+ from pyopensn.solver import DiscreteOrdinatesProblem, SteadyStateSourceSolver
+ from pyopensn.source import VolumetricSource
+ from pyopensn.xs import MultiGroupXS
+
+if __name__ == "__main__":
+ # Setup mesh
+ meshgen = FromFileMeshGenerator(
+ filename="../../../../assets/mesh/square_mesh2x2_quads_block.obj",
+ partitioner=KBAGraphPartitioner(
+ nx=2,
+ ny=2,
+ nz=1,
+ xcuts=[0.0],
+ ycuts=[0.0],
+ ),
+ )
+ grid = meshgen.Execute()
+ grid.SetOrthogonalBoundaries()
+
+ # Cross-section data
+ vol0 = RPPLogicalVolume(infx=True, infy=True, infz=True)
+ grid.SetBlockIDFromLogicalVolume(vol0, 0, True)
+ num_groups = 168
+ xs_3_170 = MultiGroupXS()
+ xs_3_170.LoadFromOpenSn("../../../../assets/xs/xs_168g.xs")
+
+ # Volumetric sources
+ strength = []
+ for g in range(num_groups):
+ strength.append(0.0)
+ mg_src1 = VolumetricSource(block_ids=[1], group_strength=strength)
+ mg_src2 = VolumetricSource(block_ids=[2], group_strength=strength)
+
+ # Boundary sources
+ bsrc = []
+ for g in range(num_groups):
+ bsrc.append(0.0)
+ bsrc[0] = 1.0
+
+ # Angular quadrature
+ pquad = GLCProductQuadrature2DXY(n_polar=2, n_azimuthal=8, scattering_order=1)
+
+ # Create solver
+ phys = DiscreteOrdinatesProblem(
+ mesh=grid,
+ num_groups=num_groups,
+ groupsets=[
+ {
+ "groups_from_to": (0, 62),
+ "angular_quadrature": pquad,
+ "angle_aggregation_num_subsets": 1,
+ "inner_linear_method": "petsc_gmres",
+ "l_abs_tol": 1.0e-6,
+ "l_max_its": 300,
+ "gmres_restart_interval": 100,
+ },
+ {
+ "groups_from_to": (63, num_groups - 1),
+ "angular_quadrature": pquad,
+ "angle_aggregation_num_subsets": 1,
+ "inner_linear_method": "petsc_gmres",
+ "l_abs_tol": 1.0e-6,
+ "l_max_its": 300,
+ "gmres_restart_interval": 100,
+ },
+ ],
+ xs_map=[{"block_ids": [0], "xs": xs_3_170}],
+ volumetric_sources=[mg_src1, mg_src2],
+ boundary_conditions=[
+ {"name": "xmin", "type": "isotropic", "group_strength": bsrc},
+ ],
+ options={"max_ags_iterations": 1},
+ sweep_type="CBC",
+ )
+ ss_solver = SteadyStateSourceSolver(problem=phys)
+ ss_solver.Initialize()
+ ss_solver.Execute()
+
+ # PPS over xmin boundary
+ pps_xmin_min = SurfacePostprocessor(
+ problem=phys,
+ value_type="min",
+ current_type="net",
+ boundaries=["xmin"]
+ )
+ pps_xmin_min.Execute()
+
+ min = pps_xmin_min.GetValue()[0]
+ assert math.isclose(min[0], -0.3341564, rel_tol=1e-6)
+ assert math.isclose(min[1], 2.09176149e-4, rel_tol=1e-6)
+
+ pps_xmin_max = SurfacePostprocessor(
+ problem=phys,
+ value_type="max",
+ current_type="net",
+ boundaries=["xmin"]
+ )
+ pps_xmin_max.Execute()
+
+ max = pps_xmin_max.GetValue()[0]
+ assert math.isclose(max[0], -0.33314547, rel_tol=1e-6)
+ assert math.isclose(max[1], 1.3370162e-3, rel_tol=1e-6)
+
+ pps_xmin_int = SurfacePostprocessor(
+ problem=phys,
+ value_type="integral",
+ current_type="net",
+ boundaries=["xmin"]
+ )
+ pps_xmin_int.Execute()
+
+ int = pps_xmin_int.GetValue()[0]
+ assert math.isclose(int[0], -10.67881732, rel_tol=1e-6)
+ assert math.isclose(int[1], 2.723323766e-2, rel_tol=1e-6)
+
+ # PPS over xmin boundary for a single groups
+ pps_xmin_g1_min = SurfacePostprocessor(
+ problem=phys,
+ value_type="min",
+ group=4,
+ current_type="net",
+ boundaries=["xmin"]
+ )
+ pps_xmin_g1_min.Execute()
+
+ min = pps_xmin_g1_min.GetValue()[0]
+ assert math.isclose(min[0], 4.49308817e-4, rel_tol=1e-6)
+
+ pps_xmin_g1_max = SurfacePostprocessor(
+ problem=phys,
+ value_type="max",
+ group=4,
+ current_type="net",
+ boundaries=["xmin"]
+ )
+ pps_xmin_g1_max.Execute()
+
+ max = pps_xmin_g1_max.GetValue()[0]
+ assert math.isclose(max[0], 6.39729658e-4, rel_tol=1e-6)
+
+ pps_xmin_g1_int = SurfacePostprocessor(
+ problem=phys,
+ value_type="integral",
+ group=4,
+ current_type="net",
+ boundaries=["xmin"]
+ )
+ pps_xmin_g1_int.Execute()
+
+ int = pps_xmin_g1_int.GetValue()[0]
+ assert math.isclose(int[0], 1.81235559e-2, rel_tol=1e-6)
+
+ # PPS over xmin boundary for a groupset 1
+ pps_xmin_gs1_min = SurfacePostprocessor(
+ problem=phys,
+ value_type="min",
+ groupset=1,
+ current_type="net",
+ boundaries=["xmin"]
+ )
+ pps_xmin_gs1_min.Execute()
+
+ min = pps_xmin_gs1_min.GetValue()[0]
+ assert math.isclose(min[0], 9.679619460480029e-05, rel_tol=1e-6)
+ assert math.isclose(min[1], 8.67902372833718e-05, rel_tol=1e-6)
+
+ pps_xmin_gs1_max = SurfacePostprocessor(
+ problem=phys,
+ value_type="max",
+ groupset=1,
+ current_type="net",
+ boundaries=["xmin"]
+ )
+ pps_xmin_gs1_max.Execute()
+
+ max = pps_xmin_gs1_max.GetValue()[0]
+ assert math.isclose(max[0], 5.97642956e-4, rel_tol=1e-6)
+ assert math.isclose(max[1], 5.38253292e-4, rel_tol=1e-6)
+
+ pps_xmin_gs1_int = SurfacePostprocessor(
+ problem=phys,
+ value_type="integral",
+ groupset=1,
+ current_type="net",
+ boundaries=["xmin"]
+ )
+ pps_xmin_gs1_int.Execute()
+
+ int = pps_xmin_gs1_int.GetValue()[0]
+ assert math.isclose(int[0], 1.32020085e-2, rel_tol=1e-6)
+ assert math.isclose(int[1], 1.18703524e-2, rel_tol=1e-6)