CBC rz sweep chunk

doc/source/userguide/discrete_ordinates_problems.rst

@@ -303,6 +303,8 @@ Practical recommendation:
 * ``AAH`` remains the default production choice and is the safer option for
   most users, particularly for problems with cyclic sweep dependencies.
 * Both ``AAH`` and ``CBC`` support time-dependent (transient) mode.
+* Both ``AAH`` and ``CBC`` are available for CPU curvilinear
+  :py:class:`pyopensn.solver.DiscreteOrdinatesCurvilinearProblem` solves.
 * Choose ``CBC`` only when the sweep graph is known to be acyclic or when
   you have verified it meets the acyclicity requirement for your specific
   problem.
@@ -524,9 +526,8 @@ curvilinear companion to the Cartesian problem class.
 
 It uses the same general construction pattern, but currently requires:
 
-* a suitable curvilinear mesh,
-* ``coord_system=2`` for cylindrical coordinates,
-* a compatible quadrature and solver setup.
+* a suitable curvilinear mesh
+* ``coord_system=2`` for cylindrical coordinates
 
 Important current limitations:
 
@@ -548,6 +549,12 @@ Example:
        sweep_type="AAH",
    )
 
+For CPU curvilinear solves, ``sweep_type`` may be either ``"AAH"`` or
+``"CBC"``. The same sweep-cycle guidance applies as for Cartesian problems:
+use ``AAH`` as the default when cyclic sweep dependencies are possible, and
+choose ``CBC`` only for problems whose sweep graph satisfies CBC's acyclicity
+requirements.
+
 Typical Construction Patterns
 =============================
 

doc/source/userguide/example_problems.rst

@@ -541,6 +541,11 @@ problem type and a compatible quadrature.
    solver.Initialize()
    solver.Execute()
 
+The curvilinear problem supports CPU ``AAH`` and ``CBC`` sweep types. Use
+``AAH`` as the default, particularly when cyclic sweep dependencies are
+possible; choose ``CBC`` only when the sweep graph satisfies CBC's acyclicity
+requirements.
+
 Example 9: Updating a Problem In Place
 ======================================
 

doc/source/userguide/faq.rst

@@ -90,6 +90,10 @@ Both ``AAH`` and ``CBC`` support time-dependent (transient) mode. Choose ``CBC``
 only if the sweep graph is known to be acyclic or has been verified to meet
 ``CBC``'s acyclicity requirements for your specific problem.
 
+For CPU curvilinear problems, both ``AAH`` and ``CBC`` are available. The same
+acyclicity requirement applies when using ``CBC``. Curvilinear GPU sweeps are
+not yet supported.
+
 Do I need to export field functions during every run?
 =====================================================
 

doc/source/userguide/solvers.rst

@@ -312,6 +312,12 @@ curvilinear geometry and currently requires:
 This problem type is experimental, and it should be treated that way in
 production workflows.
 
+CPU curvilinear solves support both ``sweep_type="AAH"`` and
+``sweep_type="CBC"``. Use ``AAH`` as the default, particularly when
+cyclic sweep dependencies are possible. Choose ``CBC`` only when the sweep graph
+is acyclic for the chosen mesh, decomposition, and quadrature. Curvilinear GPU
+sweeps are not supported.
+
 .. note::
 
    The curvilinear problem class is best used when the mesh and quadrature are
@@ -336,6 +342,7 @@ Constructor
 The constructor takes:
 
 * ``problem``: an existing :py:class:`pyopensn.solver.DiscreteOrdinatesProblem`
+  or :py:class:`pyopensn.solver.DiscreteOrdinatesCurvilinearProblem`
 
 GPU support
 -----------
@@ -688,7 +695,8 @@ Constructor
 The constructor takes:
 
 * ``problem``: an existing
-  :py:class:`pyopensn.solver.DiscreteOrdinatesProblem`
+  :py:class:`pyopensn.solver.DiscreteOrdinatesProblem` or
+  :py:class:`pyopensn.solver.DiscreteOrdinatesCurvilinearProblem`
 * ``max_iters``: maximum power iterations
 * ``k_tol``: convergence tolerance on ``k_eff``
 * ``reset_phi0``: whether to reset scalar fluxes to 1.0 before solving
@@ -770,6 +778,7 @@ Constructor
 The constructor takes:
 
 * ``problem``: an existing :py:class:`pyopensn.solver.DiscreteOrdinatesProblem`
+  or :py:class:`pyopensn.solver.DiscreteOrdinatesCurvilinearProblem`
 * nonlinear tolerances:
   ``nl_abs_tol``, ``nl_rel_tol``, ``nl_sol_tol``, ``nl_max_its``
 * linear tolerances:

modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/discrete_ordinates_curvilinear_problem.cc

@@ -3,6 +3,7 @@
 
 #include "modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/discrete_ordinates_curvilinear_problem.h"
 #include "modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/sweep_chunks/aah_sweep_chunk_rz.h"
+#include "modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/sweep_chunks/cbc_sweep_chunk_rz.h"
 #include "framework/math/spatial_discretization/finite_element/piecewise_linear/piecewise_linear_discontinuous.h"
 #include "framework/math/quadratures/angular/curvilinear_product_quadrature.h"
 #include "framework/mesh/mesh_continuum/mesh_continuum.h"
@@ -101,7 +102,7 @@ DiscreteOrdinatesCurvilinearProblem::PerformInputChecks()
       "geometries yet.");
   }
 
-  for (size_t gs = 0; gs < groupsets_.size(); ++gs)
+  for (std::size_t gs = 0; gs < groupsets_.size(); ++gs)
   {
     // angular quadrature type must be compatible with coordinate system
     const auto angular_quad_ptr = groupsets_[gs].quadrature;
@@ -209,7 +210,7 @@ DiscreteOrdinatesCurvilinearProblem::PerformInputChecks()
       if (not face.has_neighbor)
       {
         bool face_orthogonal = false;
-        for (size_t d = 0; d < unit_normal_vectors.size(); ++d)
+        for (std::size_t d = 0; d < unit_normal_vectors.size(); ++d)
         {
           const auto n_dot_e = face.normal.Dot(unit_normal_vectors[d]);
           if (std::fabs(n_dot_e) > 0.999999)
@@ -289,8 +290,7 @@ DiscreteOrdinatesCurvilinearProblem::ComputeSecondaryUnitIntegrals()
   auto ComputeCellUnitIntegrals = [&sdm, &swf](const Cell& cell)
   {
     const auto& cell_mapping = sdm.GetCellMapping(cell);
-    //    const size_t cell_num_faces = cell.faces.size();
-    const size_t cell_num_nodes = cell_mapping.GetNumNodes();
+    const std::size_t cell_num_nodes = cell_mapping.GetNumNodes();
     const auto fe_vol_data = cell_mapping.MakeVolumetricFiniteElementData();
 
     DenseMatrix<double> IntV_shapeI_shapeJ(cell_num_nodes, cell_num_nodes, 0.0);
@@ -319,7 +319,7 @@ DiscreteOrdinatesCurvilinearProblem::ComputeSecondaryUnitIntegrals()
                             {}};
   };
 
-  const size_t num_local_cells = grid_->local_cells.size();
+  const std::size_t num_local_cells = grid_->local_cells.size();
   secondary_unit_cell_matrices_.resize(num_local_cells);
 
   for (const auto& cell : grid_->local_cells)
@@ -338,9 +338,13 @@ DiscreteOrdinatesCurvilinearProblem::GetSecondaryUnitCellMatrices() const
 std::shared_ptr<SweepChunk>
 DiscreteOrdinatesCurvilinearProblem::SetSweepChunk(LBSGroupset& groupset)
 {
-  auto sweep_chunk = std::make_shared<AAHSweepChunkRZ>(*this, groupset);
+  if (sweep_type_ == "AAH")
+    return std::make_shared<AAHSweepChunkRZ>(*this, groupset);
 
-  return sweep_chunk;
+  if (sweep_type_ == "CBC")
+    return std::make_shared<CBCSweepChunkRZ>(*this, groupset);
+
+  OpenSnLogicalError("Unsupported sweep_type_ \"" + sweep_type_ + "\"");
 }
 
 } // namespace opensn