diff --git a/doc/source/userguide/discrete_ordinates_problems.rst b/doc/source/userguide/discrete_ordinates_problems.rst index 7f9b1bc99..ea8e7e4f8 100644 --- a/doc/source/userguide/discrete_ordinates_problems.rst +++ b/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 ============================= diff --git a/doc/source/userguide/example_problems.rst b/doc/source/userguide/example_problems.rst index df76747ca..b08818917 100644 --- a/doc/source/userguide/example_problems.rst +++ b/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 ====================================== diff --git a/doc/source/userguide/faq.rst b/doc/source/userguide/faq.rst index 9c49f7b82..c392e32ba 100644 --- a/doc/source/userguide/faq.rst +++ b/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? ===================================================== diff --git a/doc/source/userguide/solvers.rst b/doc/source/userguide/solvers.rst index fbcc05dcb..94c6d06a1 100644 --- a/doc/source/userguide/solvers.rst +++ b/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: diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/discrete_ordinates_curvilinear_problem.cc b/modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/discrete_ordinates_curvilinear_problem.cc index 51e665a4b..5b7ce6de8 100644 --- a/modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/discrete_ordinates_curvilinear_problem.cc +++ b/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 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 DiscreteOrdinatesCurvilinearProblem::SetSweepChunk(LBSGroupset& groupset) { - auto sweep_chunk = std::make_shared(*this, groupset); + if (sweep_type_ == "AAH") + return std::make_shared(*this, groupset); - return sweep_chunk; + if (sweep_type_ == "CBC") + return std::make_shared(*this, groupset); + + OpenSnLogicalError("Unsupported sweep_type_ \"" + sweep_type_ + "\""); } } // namespace opensn diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/sweep_chunks/aah_sweep_chunk_rz.cc b/modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/sweep_chunks/aah_sweep_chunk_rz.cc index 4b162193b..a9e75f366 100644 --- a/modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/sweep_chunks/aah_sweep_chunk_rz.cc +++ b/modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/sweep_chunks/aah_sweep_chunk_rz.cc @@ -4,10 +4,14 @@ #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_problem/sweep/fluds/aah_fluds.h" +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep_chunks/avx_sweep_chunk_utils.h" #include "modules/linear_boltzmann_solvers/lbs_problem/groupset/lbs_groupset.h" #include "framework/math/spatial_discretization/spatial_discretization.h" #include "framework/math/quadratures/angular/curvilinear_product_quadrature.h" #include "framework/mesh/mesh_continuum/mesh_continuum.h" +#include "framework/utils/error.h" +#include +#include #include namespace opensn @@ -31,7 +35,9 @@ AAHSweepChunkRZ::AAHSweepChunkRZ(DiscreteOrdinatesProblem& problem, LBSGroupset& .GetSecondaryUnitCellMatrices()), unknown_manager_(), psi_sweep_(), - normal_vector_boundary_() + normal_vector_boundary_(), + group_block_size_(ComputeGroupBlockSize(groupset.GetNumGroups())), + sweep_impl_(&AAHSweepChunkRZ::Sweep_Generic) { const auto curvilinear_product_quadrature = std::dynamic_pointer_cast(groupset_.quadrature); @@ -41,8 +47,8 @@ AAHSweepChunkRZ::AAHSweepChunkRZ(DiscreteOrdinatesProblem& problem, LBSGroupset& "invalid angular quadrature"); // configure unknown manager for quantities that depend on polar level - const size_t dir_map_size = curvilinear_product_quadrature->GetDirectionMap().size(); - for (size_t m = 0; m < dir_map_size; ++m) + const std::size_t dir_map_size = curvilinear_product_quadrature->GetDirectionMap().size(); + for (std::size_t m = 0; m < dir_map_size; ++m) unknown_manager_.AddUnknown(UnknownType::VECTOR_N, groupset_.GetNumGroups()); // allocate storage for sweeping dependency @@ -58,10 +64,46 @@ AAHSweepChunkRZ::AAHSweepChunkRZ(DiscreteOrdinatesProblem& problem, LBSGroupset& const auto d = (grid_->GetDimension() == 1) ? 2 : 0; normal_vector_boundary_ = Vector3(0.0, 0.0, 0.0); normal_vector_boundary_(d) = 1; + + if (min_num_cell_dofs_ == max_num_cell_dofs_) + { + switch (min_num_cell_dofs_) + { + case 2: + sweep_impl_ = &AAHSweepChunkRZ::Sweep_FixedN<2>; + break; + case 3: + sweep_impl_ = &AAHSweepChunkRZ::Sweep_FixedN<3>; + break; + case 4: + sweep_impl_ = &AAHSweepChunkRZ::Sweep_FixedN<4>; + break; + case 5: + sweep_impl_ = &AAHSweepChunkRZ::Sweep_FixedN<5>; + break; + case 6: + sweep_impl_ = &AAHSweepChunkRZ::Sweep_FixedN<6>; + break; + case 7: + sweep_impl_ = &AAHSweepChunkRZ::Sweep_FixedN<7>; + break; + case 8: + sweep_impl_ = &AAHSweepChunkRZ::Sweep_FixedN<8>; + break; + default: + break; + } + } } void AAHSweepChunkRZ::Sweep(AngleSet& angle_set) +{ + (this->*sweep_impl_)(angle_set); +} + +void +AAHSweepChunkRZ::Sweep_Generic(AngleSet& angle_set) { auto gs_size = groupset_.GetNumGroups(); auto gs_gi = groupset_.first_group; @@ -84,8 +126,8 @@ AAHSweepChunkRZ::Sweep(AngleSet& angle_set) // Loop over each cell const auto& spds = angle_set.GetSPDS(); const auto& spls = spds.GetLocalSubgrid(); - const size_t num_spls = spls.size(); - for (size_t spls_index = 0; spls_index < num_spls; ++spls_index) + const std::size_t num_spls = spls.size(); + for (std::size_t spls_index = 0; spls_index < num_spls; ++spls_index) { auto cell_local_id = spls[spls_index]; const auto& cell = grid_->local_cells[cell_local_id]; @@ -110,7 +152,7 @@ AAHSweepChunkRZ::Sweep(AngleSet& angle_set) const auto ni_deploc_face_counter = deploc_face_counter; const auto ni_preloc_face_counter = preloc_face_counter; const std::vector& as_angle_indices = angle_set.GetAngleIndices(); - for (size_t as_ss_idx = 0; as_ss_idx < as_angle_indices.size(); ++as_ss_idx) + for (std::size_t as_ss_idx = 0; as_ss_idx < as_angle_indices.size(); ++as_ss_idx) { auto direction_num = as_angle_indices[as_ss_idx]; auto omega = groupset_.quadrature->GetOmega(direction_num); @@ -126,31 +168,31 @@ AAHSweepChunkRZ::Sweep(AngleSet& angle_set) preloc_face_counter = ni_preloc_face_counter; // Reset right-hand side - for (size_t gsg = 0; gsg < gs_size; ++gsg) + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) b[gsg] = Vector(cell_num_nodes, 0.0); - for (size_t i = 0; i < cell_num_nodes; ++i) + for (std::size_t i = 0; i < cell_num_nodes; ++i) { - for (size_t j = 0; j < cell_num_nodes; ++j) + for (std::size_t j = 0; j < cell_num_nodes; ++j) { const auto jr = discretization_.MapDOFLocal(cell, j, unknown_manager_, polar_level, gs_gi); - for (size_t gsg = 0; gsg < gs_size; ++gsg) + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) b[gsg](i) += fac_streaming_operator * Maux(i, j) * psi_sweep_[jr + gsg]; } } - for (size_t i = 0; i < cell_num_nodes; ++i) - for (size_t j = 0; j < cell_num_nodes; ++j) + for (std::size_t i = 0; i < cell_num_nodes; ++i) + for (std::size_t j = 0; j < cell_num_nodes; ++j) Amat(i, j) = omega.Dot(G(i, j)) + fac_streaming_operator * Maux(i, j); // Update face orientations - for (size_t f = 0; f < cell_num_faces; ++f) + for (std::size_t f = 0; f < cell_num_faces; ++f) face_mu_values[f] = omega.Dot(cell.faces[f].normal); // Surface integrals int in_face_counter = -1; - for (size_t f = 0; f < cell_num_faces; ++f) + for (std::size_t f = 0; f < cell_num_faces; ++f) { if (face_orientations[f] != FaceOrientation::INCOMING) continue; @@ -165,12 +207,12 @@ AAHSweepChunkRZ::Sweep(AngleSet& angle_set) ++preloc_face_counter; // IntSf_mu_psi_Mij_dA - const size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f); - for (size_t fi = 0; fi < num_face_nodes; ++fi) + const std::size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f); + for (std::size_t fi = 0; fi < num_face_nodes; ++fi) { const int i = cell_mapping.MapFaceNode(f, fi); - for (size_t fj = 0; fj < num_face_nodes; ++fj) + for (std::size_t fj = 0; fj < num_face_nodes; ++fj) { const int j = cell_mapping.MapFaceNode(f, fj); @@ -210,24 +252,24 @@ AAHSweepChunkRZ::Sweep(AngleSet& angle_set) if (not psi) continue; - for (size_t gsg = 0; gsg < gs_size; ++gsg) + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) b[gsg](i) += psi[gsg] * mu_Nij; } // for face node j } // for face node i } // for f const auto row_offset = - static_cast(direction_num) * static_cast(num_moments_); + static_cast(direction_num) * static_cast(num_moments_); const double* m2d_row = m2d_op.data() + row_offset; const double* d2m_row = d2m_op.data() + row_offset; // Looping over groups, assembling mass terms - for (size_t gsg = 0; gsg < gs_size; ++gsg) + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) { double sigma_tg = sigma_t[gs_gi + gsg]; // Contribute source moments q = M_n^T * q_moms - for (size_t i = 0; i < cell_num_nodes; ++i) + for (std::size_t i = 0; i < cell_num_nodes; ++i) { double temp_src = 0.0; for (unsigned int m = 0; m < num_moments_; ++m) @@ -241,10 +283,10 @@ AAHSweepChunkRZ::Sweep(AngleSet& angle_set) // Mass matrix and source // Atemp = Amat + sigma_tgr * M // b += M * q - for (size_t i = 0; i < cell_num_nodes; ++i) + for (std::size_t i = 0; i < cell_num_nodes; ++i) { double temp = 0.0; - for (size_t j = 0; j < cell_num_nodes; ++j) + for (std::size_t j = 0; j < cell_num_nodes; ++j) { const double Mij = M(i, j); Atemp(i, j) = Amat(i, j) + Mij * sigma_tg; @@ -261,10 +303,10 @@ AAHSweepChunkRZ::Sweep(AngleSet& angle_set) for (unsigned int m = 0; m < num_moments_; ++m) { const double wn_d2m = d2m_row[m]; - for (size_t i = 0; i < cell_num_nodes; ++i) + for (std::size_t i = 0; i < cell_num_nodes; ++i) { const auto ir = cell_transport_view.MapDOF(i, m, gs_gi); - for (size_t gsg = 0; gsg < gs_size; ++gsg) + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) destination_phi_[ir + gsg] += wn_d2m * b[gsg](i); } } @@ -275,11 +317,11 @@ AAHSweepChunkRZ::Sweep(AngleSet& angle_set) double* cell_psi_data = &destination_psi_[discretization_.MapDOFLocal(cell, 0, groupset_.psi_uk_man_, 0, 0)]; - for (size_t i = 0; i < cell_num_nodes; ++i) + for (std::size_t i = 0; i < cell_num_nodes; ++i) { - const size_t imap = + const std::size_t imap = i * groupset_angle_group_stride_ + direction_num * groupset_group_stride_; - for (size_t gsg = 0; gsg < gs_size; ++gsg) + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) cell_psi_data[imap + gsg] = b[gsg](i); } } @@ -287,7 +329,7 @@ AAHSweepChunkRZ::Sweep(AngleSet& angle_set) // For outgoing, non-boundary faces, copy angular flux to fluds and // accumulate outflow int out_face_counter = -1; - for (size_t f = 0; f < cell_num_faces; ++f) + for (std::size_t f = 0; f < cell_num_faces; ++f) { if (face_orientations[f] != FaceOrientation::OUTGOING) continue; @@ -303,14 +345,14 @@ AAHSweepChunkRZ::Sweep(AngleSet& angle_set) if (not is_boundary_face and not is_local_face) ++deploc_face_counter; - const size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f); - for (size_t fi = 0; fi < num_face_nodes; ++fi) + const std::size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f); + for (std::size_t fi = 0; fi < num_face_nodes; ++fi) { const int i = cell_mapping.MapFaceNode(f, fi); if (is_boundary_face) { - for (size_t gsg = 0; gsg < gs_size; ++gsg) + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) cell_outflow_view.Add( f, gs_gi + gsg, wt * face_mu_values[f] * b[gsg](i) * IntF_shapeI(i)); } @@ -327,7 +369,7 @@ AAHSweepChunkRZ::Sweep(AngleSet& angle_set) if (not is_boundary_face or is_reflecting_boundary_face) { - for (size_t gsg = 0; gsg < gs_size; ++gsg) + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) psi[gsg] = b[gsg](i); } } // for fi @@ -337,14 +379,362 @@ AAHSweepChunkRZ::Sweep(AngleSet& angle_set) // interval) const auto f0 = 1 / fac_diamond_difference; const auto f1 = f0 - 1; - for (size_t i = 0; i < cell_num_nodes; ++i) + for (std::size_t i = 0; i < cell_num_nodes; ++i) { const auto ir = discretization_.MapDOFLocal(cell, i, unknown_manager_, polar_level, gs_gi); - for (size_t gsg = 0; gsg < gs_size; ++gsg) + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) psi_sweep_[ir + gsg] = f0 * b[gsg](i) - f1 * psi_sweep_[ir + gsg]; } } // for angleset/subset } // for cell } +template +void +AAHSweepChunkRZ::Sweep_FixedN(AngleSet& angle_set) +{ + static_assert(NumNodes >= 2 and NumNodes <= 8); + + const auto gs_size = groupset_.GetNumGroups(); + const auto gs_gi = groupset_.first_group; + + int deploc_face_counter = -1; + int preloc_face_counter = -1; + + auto& fluds = dynamic_cast(angle_set.GetFLUDS()); + const auto& m2d_op = groupset_.quadrature->GetMomentToDiscreteOperator(); + const auto& d2m_op = groupset_.quadrature->GetDiscreteToMomentOperator(); + + constexpr std::size_t matrix_size = static_cast(NumNodes) * NumNodes; + constexpr auto idx = [](std::size_t i, std::size_t j) -> std::size_t { return i * NumNodes + j; }; + + std::vector b(static_cast(gs_size) * NumNodes, 0.0); + std::vector sigma_block; + sigma_block.reserve(group_block_size_); + std::array source{}; + std::vector> moment_dof_map(num_moments_); + + const auto curvilinear_product_quadrature = + std::dynamic_pointer_cast(groupset_.quadrature); + OpenSnLogicalErrorIf(curvilinear_product_quadrature == nullptr, + "AAHSweepChunkRZ requires a curvilinear product quadrature."); + + const auto& spds = angle_set.GetSPDS(); + const auto& spls = spds.GetLocalSubgrid(); + const std::size_t num_spls = spls.size(); + for (std::size_t spls_index = 0; spls_index < num_spls; ++spls_index) + { + const auto cell_local_id = spls[spls_index]; + const auto& cell = grid_->local_cells[cell_local_id]; + const auto& cell_mapping = discretization_.GetCellMapping(cell); + const auto& cell_transport_view = cell_transport_views_[cell_local_id]; + auto& cell_outflow_view = cell_outflow_views_[cell_local_id]; + const auto cell_num_faces = cell.faces.size(); + + OpenSnInvalidArgumentIf(cell_mapping.GetNumNodes() != static_cast(NumNodes), + "AAHSweepChunkRZ::Sweep_FixedN invoked for an incompatible cell " + "topology."); + + const auto& face_orientations = spds.GetCellFaceOrientations()[cell_local_id]; + std::vector face_mu_values(cell_num_faces); + + const auto& sigma_t = xs_.at(cell.block_id)->GetSigmaTotal(); + + const auto& G = unit_cell_matrices_[cell_local_id].intV_shapeI_gradshapeJ; + const auto& M = unit_cell_matrices_[cell_local_id].intV_shapeI_shapeJ; + const auto& M_surf = unit_cell_matrices_[cell_local_id].intS_shapeI_shapeJ; + const auto& Maux = secondary_unit_cell_matrices_[cell_local_id].intV_shapeI_shapeJ; + + std::array mass_matrix{}; + std::array aux_matrix{}; + PRAGMA_UNROLL + for (std::size_t i = 0; i < NumNodes; ++i) + { + PRAGMA_UNROLL + for (std::size_t j = 0; j < NumNodes; ++j) + { + mass_matrix[idx(i, j)] = M(i, j); + aux_matrix[idx(i, j)] = Maux(i, j); + } + } + + for (unsigned int m = 0; m < num_moments_; ++m) + { + PRAGMA_UNROLL + for (std::size_t i = 0; i < NumNodes; ++i) + moment_dof_map[m][i] = cell_transport_view.MapDOF(i, m, gs_gi); + } + + const auto ni_deploc_face_counter = deploc_face_counter; + const auto ni_preloc_face_counter = preloc_face_counter; + const auto& as_angle_indices = angle_set.GetAngleIndices(); + for (std::size_t as_ss_idx = 0; as_ss_idx < as_angle_indices.size(); ++as_ss_idx) + { + const auto direction_num = as_angle_indices[as_ss_idx]; + const auto& omega = groupset_.quadrature->GetOmega(direction_num); + const auto wt = groupset_.quadrature->GetWeight(direction_num); + + const auto polar_level = map_polar_level_[direction_num]; + const auto fac_diamond_difference = + curvilinear_product_quadrature->GetDiamondDifferenceFactor()[direction_num]; + const auto fac_streaming_operator = + curvilinear_product_quadrature->GetStreamingOperatorFactor()[direction_num]; + + deploc_face_counter = ni_deploc_face_counter; + preloc_face_counter = ni_preloc_face_counter; + + std::fill(b.begin(), b.end(), 0.0); + + std::array Amat{}; + PRAGMA_UNROLL + for (std::size_t i = 0; i < NumNodes; ++i) + { + PRAGMA_UNROLL + for (std::size_t j = 0; j < NumNodes; ++j) + { + const auto jr = + discretization_.MapDOFLocal(cell, j, unknown_manager_, polar_level, gs_gi); + const double streaming = fac_streaming_operator * aux_matrix[idx(i, j)]; + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) + b[gsg * NumNodes + i] += streaming * psi_sweep_[jr + gsg]; + + Amat[idx(i, j)] = omega.Dot(G(i, j)) + streaming; + } + } + + for (std::size_t f = 0; f < cell_num_faces; ++f) + face_mu_values[f] = omega.Dot(cell.faces[f].normal); + + int in_face_counter = -1; + for (std::size_t f = 0; f < cell_num_faces; ++f) + { + if (face_orientations[f] != FaceOrientation::INCOMING) + continue; + + const auto& cell_face = cell.faces[f]; + const bool is_local_face = cell_transport_view.IsFaceLocal(f); + const bool is_boundary_face = not cell_face.has_neighbor; + + if (is_local_face) + ++in_face_counter; + else if (not is_boundary_face) + ++preloc_face_counter; + + const std::size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f); + for (std::size_t fi = 0; fi < num_face_nodes; ++fi) + { + const int i = cell_mapping.MapFaceNode(f, fi); + for (std::size_t fj = 0; fj < num_face_nodes; ++fj) + { + const int j = cell_mapping.MapFaceNode(f, fj); + const double mu_Nij = -face_mu_values[f] * M_surf[f](i, j); + Amat[idx(i, j)] += mu_Nij; + + const double* psi = nullptr; + if (is_local_face) + psi = fluds.UpwindPsi(spls_index, in_face_counter, fj, 0, as_ss_idx); + else if (not is_boundary_face) + psi = fluds.NLUpwindPsi(preloc_face_counter, fj, 0, as_ss_idx); + else + { + const bool incident_on_symmetric_boundary = + (cell_face.normal.Dot(normal_vector_boundary_) < -0.999999); + if (not incident_on_symmetric_boundary) + { + psi = angle_set.PsiBoundary(cell_face.neighbor_id, + direction_num, + cell_local_id, + f, + fj, + gs_gi, + IsSurfaceSourceActive()); + } + } + + if (psi != nullptr) + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) + b[gsg * NumNodes + i] += psi[gsg] * mu_Nij; + } + } + } + + const auto row_offset = + static_cast(direction_num) * static_cast(num_moments_); + const double* __restrict m2d_row = m2d_op.data() + row_offset; + const double* __restrict d2m_row = d2m_op.data() + row_offset; + + for (unsigned int g0 = 0; g0 < gs_size; g0 += group_block_size_) + { + const auto g1 = std::min(g0 + group_block_size_, gs_size); + const auto block_len = g1 - g0; + sigma_block.resize(block_len); + + for (unsigned int gsg = g0; gsg < g1; ++gsg) + { + const auto rel = gsg - g0; + sigma_block[rel] = sigma_t[gs_gi + gsg]; + + source.fill(0.0); + for (unsigned int m = 0; m < num_moments_; ++m) + { + const double w = m2d_row[m]; + PRAGMA_UNROLL + for (std::size_t i = 0; i < NumNodes; ++i) + source[i] += w * source_moments_[moment_dof_map[m][i] + gsg]; + } + + auto* __restrict bg = &b[static_cast(gsg) * NumNodes]; + PRAGMA_UNROLL + for (std::size_t i = 0; i < NumNodes; ++i) + { + double temp = 0.0; + PRAGMA_UNROLL + for (std::size_t j = 0; j < NumNodes; ++j) + temp += mass_matrix[idx(i, j)] * source[j]; + bg[i] += temp; + } + } + + std::size_t k = 0; +#if __AVX512F__ + for (; k + simd_width <= block_len; k += simd_width) + detail::SimdBatchSolve( + Amat.data(), mass_matrix.data(), &sigma_block[k], &b[(g0 + k) * NumNodes]); +#elif __AVX2__ + for (; k + simd_width <= block_len; k += simd_width) + detail::SimdBatchSolve( + Amat.data(), mass_matrix.data(), &sigma_block[k], &b[(g0 + k) * NumNodes]); +#endif + + for (; k < block_len; ++k) + { + const std::size_t gsg = g0 + k; + const double sigma_tg = sigma_block[k]; + std::array A{}; + PRAGMA_UNROLL + for (std::size_t i = 0; i < NumNodes; ++i) + { + PRAGMA_UNROLL + for (std::size_t j = 0; j < NumNodes; ++j) + A[idx(i, j)] = Amat[idx(i, j)] + sigma_tg * mass_matrix[idx(i, j)]; + } + + auto* __restrict bg = &b[gsg * NumNodes]; + for (std::size_t pivot = 0; pivot < NumNodes; ++pivot) + { + const double inv = 1.0 / A[idx(pivot, pivot)]; + for (std::size_t row = pivot + 1; row < NumNodes; ++row) + { + const double factor = A[idx(row, pivot)] * inv; + bg[row] -= factor * bg[pivot]; + PRAGMA_UNROLL + for (std::size_t col = pivot + 1; col < NumNodes; ++col) + A[idx(row, col)] -= factor * A[idx(pivot, col)]; + } + } + + for (std::size_t pivot = NumNodes; pivot-- > 0;) + { + PRAGMA_UNROLL + for (std::size_t col = pivot + 1; col < NumNodes; ++col) + bg[pivot] -= A[idx(pivot, col)] * bg[col]; + bg[pivot] /= A[idx(pivot, pivot)]; + } + } + + for (unsigned int gsg = g0; gsg < g1; ++gsg) + { + const auto* __restrict bg = &b[static_cast(gsg) * NumNodes]; + for (unsigned int m = 0; m < num_moments_; ++m) + { + const double wn_d2m = d2m_row[m]; + PRAGMA_UNROLL + for (std::size_t i = 0; i < NumNodes; ++i) + destination_phi_[moment_dof_map[m][i] + gsg] += wn_d2m * bg[i]; + } + } + } + + if (SaveAngularFluxEnabled()) + { + double* cell_psi_data = + &destination_psi_[discretization_.MapDOFLocal(cell, 0, groupset_.psi_uk_man_, 0, 0)]; + PRAGMA_UNROLL + for (std::size_t i = 0; i < NumNodes; ++i) + { + const std::size_t imap = + i * groupset_angle_group_stride_ + direction_num * groupset_group_stride_; + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) + cell_psi_data[imap + gsg] = b[gsg * NumNodes + i]; + } + } + + int out_face_counter = -1; + for (std::size_t f = 0; f < cell_num_faces; ++f) + { + if (face_orientations[f] != FaceOrientation::OUTGOING) + continue; + + ++out_face_counter; + const auto& face = cell.faces[f]; + const bool is_local_face = cell_transport_view.IsFaceLocal(f); + const bool is_boundary_face = not face.has_neighbor; + const bool is_reflecting_boundary_face = + (is_boundary_face and angle_set.GetBoundaries()[face.neighbor_id]->IsReflecting()); + const auto& IntF_shapeI = unit_cell_matrices_[cell_local_id].intS_shapeI[f]; + + if (not is_boundary_face and not is_local_face) + ++deploc_face_counter; + + const std::size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f); + for (std::size_t fi = 0; fi < num_face_nodes; ++fi) + { + const int i = cell_mapping.MapFaceNode(f, fi); + + if (is_boundary_face) + { + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) + cell_outflow_view.Add( + f, gs_gi + gsg, wt * face_mu_values[f] * b[gsg * NumNodes + i] * IntF_shapeI(i)); + } + + double* psi = nullptr; + if (is_local_face) + psi = fluds.OutgoingPsi(spls_index, out_face_counter, fi, as_ss_idx); + else if (not is_boundary_face) + psi = fluds.NLOutgoingPsi(deploc_face_counter, fi, as_ss_idx); + else if (is_reflecting_boundary_face) + psi = angle_set.PsiReflected(face.neighbor_id, direction_num, cell_local_id, f, fi); + else + continue; + + if (not is_boundary_face or is_reflecting_boundary_face) + { + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) + psi[gsg] = b[gsg * NumNodes + i]; + } + } + } + + const auto f0 = 1.0 / fac_diamond_difference; + const auto f1 = f0 - 1.0; + PRAGMA_UNROLL + for (std::size_t i = 0; i < NumNodes; ++i) + { + const auto ir = discretization_.MapDOFLocal(cell, i, unknown_manager_, polar_level, gs_gi); + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) + psi_sweep_[ir + gsg] = f0 * b[gsg * NumNodes + i] - f1 * psi_sweep_[ir + gsg]; + } + } + } +} + +template void AAHSweepChunkRZ::Sweep_FixedN<2>(AngleSet&); +template void AAHSweepChunkRZ::Sweep_FixedN<3>(AngleSet&); +template void AAHSweepChunkRZ::Sweep_FixedN<4>(AngleSet&); +template void AAHSweepChunkRZ::Sweep_FixedN<5>(AngleSet&); +template void AAHSweepChunkRZ::Sweep_FixedN<6>(AngleSet&); +template void AAHSweepChunkRZ::Sweep_FixedN<7>(AngleSet&); +template void AAHSweepChunkRZ::Sweep_FixedN<8>(AngleSet&); + } // namespace opensn diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/sweep_chunks/aah_sweep_chunk_rz.h b/modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/sweep_chunks/aah_sweep_chunk_rz.h index 9c397cff0..0b8f253aa 100644 --- a/modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/sweep_chunks/aah_sweep_chunk_rz.h +++ b/modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/sweep_chunks/aah_sweep_chunk_rz.h @@ -21,6 +21,13 @@ class AAHSweepChunkRZ : public SweepChunk void Sweep(AngleSet& angle_set) override; private: + using SweepFunc = void (AAHSweepChunkRZ::*)(AngleSet&); + + void Sweep_Generic(AngleSet& angle_set); + + template + void Sweep_FixedN(AngleSet& angle_set); + /// Secondary spatial discretization cell matrices const std::vector& secondary_unit_cell_matrices_; /// Unknown manager. @@ -31,6 +38,10 @@ class AAHSweepChunkRZ : public SweepChunk std::map map_polar_level_; /// Normal vector to determine symmetric boundary condition. Vector3 normal_vector_boundary_; + /// Number of groups solved in one block. + unsigned int group_block_size_; + /// Selected sweep implementation. + SweepFunc sweep_impl_; }; } // namespace opensn diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/sweep_chunks/cbc_sweep_chunk_rz.cc b/modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/sweep_chunks/cbc_sweep_chunk_rz.cc new file mode 100644 index 000000000..f89b29a5f --- /dev/null +++ b/modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/sweep_chunks/cbc_sweep_chunk_rz.cc @@ -0,0 +1,742 @@ +// SPDX-FileCopyrightText: 2026 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#include "modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/sweep_chunks/cbc_sweep_chunk_rz.h" +#include "modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/discrete_ordinates_curvilinear_problem.h" +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/fluds/cbc_fluds.h" +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep_chunks/avx_sweep_chunk_utils.h" +#include "modules/linear_boltzmann_solvers/lbs_problem/groupset/lbs_groupset.h" +#include "framework/math/quadratures/angular/curvilinear_product_quadrature.h" +#include "framework/math/spatial_discretization/spatial_discretization.h" +#include "framework/mesh/mesh_continuum/mesh_continuum.h" +#include "framework/utils/error.h" +#include "caliper/cali.h" +#include +#include +#include + +namespace opensn +{ +namespace +{ + +const CurvilinearProductQuadrature& +RequireCurvilinearProductQuadrature(const LBSGroupset& groupset) +{ + const auto* quadrature = + dynamic_cast(groupset.quadrature.get()); + if (quadrature == nullptr) + throw std::invalid_argument("CBCSweepChunkRZ requires a curvilinear product quadrature."); + return *quadrature; +} + +} // namespace + +CBCSweepChunkRZ::CBCSweepChunkRZ(DiscreteOrdinatesProblem& problem, LBSGroupset& groupset) + : SweepChunk(problem.GetPhiNewLocal(), + problem.GetPsiNewLocal()[groupset.id], + problem.GetGrid(), + problem.GetSpatialDiscretization(), + problem.GetUnitCellMatrices(), + problem.GetCellTransportViews(), + problem.GetCellOutflowViews(), + problem.GetQMomentsLocal(), + groupset, + problem.GetBlockID2XSMap(), + problem.GetNumMoments(), + problem.GetMaxCellDOFCount(), + problem.GetMinCellDOFCount()), + curvilinear_quadrature_(RequireCurvilinearProductQuadrature(groupset)), + secondary_unit_cell_matrices_(dynamic_cast(problem) + .GetSecondaryUnitCellMatrices()), + unknown_manager_(), + psi_sweep_(), + direction_polar_level_(groupset.quadrature->GetNumAngles(), 0), + normal_vector_boundary_(), + group_block_size_(ComputeGroupBlockSize(groupset.GetNumGroups())), + sweep_impl_(&CBCSweepChunkRZ::Sweep_Generic) +{ + const auto& direction_map = curvilinear_quadrature_.GetDirectionMap(); + const auto num_polar_levels = direction_map.size(); + for (std::size_t m = 0; m < num_polar_levels; ++m) + unknown_manager_.AddUnknown(UnknownType::VECTOR_N, groupset_.GetNumGroups()); + + psi_sweep_.assign(discretization_.GetNumLocalDOFs(unknown_manager_), 0.0); + + for (const auto& [polar_level, direction_indices] : direction_map) + { + for (const auto direction_index : direction_indices) + { + OpenSnLogicalErrorIf(direction_index >= direction_polar_level_.size(), + "CBCSweepChunkRZ received an invalid quadrature direction index."); + direction_polar_level_[direction_index] = polar_level; + } + } + + const auto d = (grid_->GetDimension() == 1) ? 2 : 0; + normal_vector_boundary_ = Vector3(0.0, 0.0, 0.0); + normal_vector_boundary_(d) = 1.0; + + if (min_num_cell_dofs_ == max_num_cell_dofs_) + { + switch (min_num_cell_dofs_) + { + case 2: + sweep_impl_ = &CBCSweepChunkRZ::Sweep_FixedN<2>; + break; + case 3: + sweep_impl_ = &CBCSweepChunkRZ::Sweep_FixedN<3>; + break; + case 4: + sweep_impl_ = &CBCSweepChunkRZ::Sweep_FixedN<4>; + break; + case 5: + sweep_impl_ = &CBCSweepChunkRZ::Sweep_FixedN<5>; + break; + case 6: + sweep_impl_ = &CBCSweepChunkRZ::Sweep_FixedN<6>; + break; + case 7: + sweep_impl_ = &CBCSweepChunkRZ::Sweep_FixedN<7>; + break; + case 8: + sweep_impl_ = &CBCSweepChunkRZ::Sweep_FixedN<8>; + break; + default: + break; + } + } +} + +void +CBCSweepChunkRZ::SetAngleSet(AngleSet& angle_set) +{ + CALI_CXX_MARK_SCOPE("CBCSweepChunkRZ::SetAngleSet"); + + fluds_ = &dynamic_cast(angle_set.GetFLUDS()); + async_comm_ = &dynamic_cast(*angle_set.GetCommunicator()); +} + +void +CBCSweepChunkRZ::Sweep(AngleSet& angle_set) +{ + (this->*sweep_impl_)(angle_set); +} + +void +CBCSweepChunkRZ::Sweep_Generic(AngleSet& angle_set) +{ + CALI_CXX_MARK_SCOPE("CBCSweepChunkRZ::Sweep"); + + const auto& groupset = groupset_; + const auto gs_size = groupset.GetNumGroups(); + const auto gs_gi = groupset.first_group; + const auto num_angles_in_as = angle_set.GetNumAngles(); + const auto group_angle_stride = gs_size * num_angles_in_as; + const auto& cell = *cell_; + const auto cell_local_id = cell.local_id; + const auto& cell_mapping = discretization_.GetCellMapping(cell); + const auto& cell_transport_view = cell_transport_views_[cell_local_id]; + auto& cell_outflow_view = cell_outflow_views_[cell_local_id]; + const std::size_t cell_num_faces = cell.faces.size(); + const std::size_t cell_num_nodes = cell_mapping.GetNumNodes(); + const auto& unit_mats = unit_cell_matrices_[cell_local_id]; + auto& fluds = *fluds_; + const auto& common_data = fluds.GetCommonData(); + const auto& G = unit_mats.intV_shapeI_gradshapeJ; + const auto& M = unit_mats.intV_shapeI_shapeJ; + const auto& M_surf = unit_mats.intS_shapeI_shapeJ; + const auto& IntS_shapeI = unit_mats.intS_shapeI; + const auto& Maux = secondary_unit_cell_matrices_[cell_local_id].intV_shapeI_shapeJ; + const auto& m2d_op = groupset.quadrature->GetMomentToDiscreteOperator(); + const auto& d2m_op = groupset.quadrature->GetDiscreteToMomentOperator(); + const auto& diamond_difference_factor = curvilinear_quadrature_.GetDiamondDifferenceFactor(); + const auto& streaming_operator_factor = curvilinear_quadrature_.GetStreamingOperatorFactor(); + const auto& face_orientations = angle_set.GetSPDS().GetCellFaceOrientations()[cell_local_id]; + const auto& sigma_t = xs_.at(cell.block_id)->GetSigmaTotal(); + + DenseMatrix Amat(max_num_cell_dofs_, max_num_cell_dofs_); + DenseMatrix Atemp(max_num_cell_dofs_, max_num_cell_dofs_); + std::vector> b(gs_size, Vector(max_num_cell_dofs_)); + std::vector source(max_num_cell_dofs_); + + auto& face_mu_values = workspace_.face_mu_values; + face_mu_values.assign(cell_num_faces, 0.0); + PrepareNonlocalOutgoingPsi(workspace_, + fluds, + cell, + cell_local_id, + cell_mapping, + cell_transport_view, + group_angle_stride, + face_orientations); + + const auto& as_angle_indices = angle_set.GetAngleIndices(); + for (std::size_t as_ss_idx = 0; as_ss_idx < num_angles_in_as; ++as_ss_idx) + { + const auto direction_num = as_angle_indices[as_ss_idx]; + const auto& omega = groupset.quadrature->GetOmega(direction_num); + const auto wt = groupset.quadrature->GetWeight(direction_num); + + const auto polar_level = direction_polar_level_[direction_num]; + const auto fac_diamond_difference = diamond_difference_factor[direction_num]; + const auto fac_streaming_operator = streaming_operator_factor[direction_num]; + + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) + for (std::size_t i = 0; i < cell_num_nodes; ++i) + b[gsg](i) = 0.0; + + for (std::size_t i = 0; i < cell_num_nodes; ++i) + { + for (std::size_t j = 0; j < cell_num_nodes; ++j) + { + const auto jr = discretization_.MapDOFLocal(cell, j, unknown_manager_, polar_level, gs_gi); + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) + b[gsg](i) += fac_streaming_operator * Maux(i, j) * psi_sweep_[jr + gsg]; + } + } + + for (std::size_t i = 0; i < cell_num_nodes; ++i) + for (std::size_t j = 0; j < cell_num_nodes; ++j) + Amat(i, j) = omega.Dot(G(i, j)) + fac_streaming_operator * Maux(i, j); + + for (std::size_t f = 0; f < cell_num_faces; ++f) + face_mu_values[f] = omega.Dot(cell.faces[f].normal); + + for (std::size_t f = 0; f < cell_num_faces; ++f) + { + if (face_orientations[f] != FaceOrientation::INCOMING) + continue; + + const auto& face = cell.faces[f]; + const bool is_boundary_face = not face.has_neighbor; + const bool is_local_face = not is_boundary_face and cell_transport_view.IsFaceLocal(f); + const auto* face_nodal_mapping = + is_boundary_face + ? nullptr + : &common_data.GetFaceNodalMapping(cell_local_id, static_cast(f)); + const auto incoming_nonlocal_slot = + (is_boundary_face or is_local_face) + ? CBC_FLUDSCommonData::INVALID_FACE_SLOT + : common_data.IncomingFaceSlot(cell_local_id, static_cast(f)); + + const std::size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f); + for (std::size_t fi = 0; fi < num_face_nodes; ++fi) + { + const int i = cell_mapping.MapFaceNode(f, fi); + + for (std::size_t fj = 0; fj < num_face_nodes; ++fj) + { + const int j = cell_mapping.MapFaceNode(f, fj); + const double mu_Nij = -face_mu_values[f] * M_surf[f](i, j); + Amat(i, j) += mu_Nij; + + const double* psi = nullptr; + if (is_local_face) + { + psi = fluds.UpwindPsi(*cell_transport_view.FaceNeighbor(f), + face_nodal_mapping->cell_node_mapping_[fj], + as_ss_idx); + } + else if (not is_boundary_face) + { + psi = fluds.NLUpwindPsi( + incoming_nonlocal_slot, face_nodal_mapping->face_node_mapping_[fj], as_ss_idx); + } + else + { + const bool incident_on_symmetric_boundary = + (face.normal.Dot(normal_vector_boundary_) < -0.999999); + if (not incident_on_symmetric_boundary) + { + psi = angle_set.PsiBoundary(face.neighbor_id, + direction_num, + cell_local_id, + static_cast(f), + static_cast(fj), + gs_gi, + IsSurfaceSourceActive()); + } + } + + if (psi != nullptr) + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) + b[gsg](i) += psi[gsg] * mu_Nij; + } + } + } + + const auto dir_moment_offset = + static_cast(direction_num) * static_cast(num_moments_); + const double* m2d_row = m2d_op.data() + dir_moment_offset; + const double* d2m_row = d2m_op.data() + dir_moment_offset; + + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) + { + const double sigma_tg = sigma_t[gs_gi + gsg]; + + for (std::size_t i = 0; i < cell_num_nodes; ++i) + { + double temp_src = 0.0; + for (unsigned int m = 0; m < num_moments_; ++m) + { + const auto ir = cell_transport_view.MapDOF(i, m, gs_gi + gsg); + temp_src += m2d_row[m] * source_moments_[ir]; + } + source[i] = temp_src; + } + + for (std::size_t i = 0; i < cell_num_nodes; ++i) + { + double temp = 0.0; + for (std::size_t j = 0; j < cell_num_nodes; ++j) + { + const double Mij = M(i, j); + Atemp(i, j) = Amat(i, j) + Mij * sigma_tg; + temp += Mij * source[j]; + } + b[gsg](i) += temp; + } + + GaussElimination(Atemp, b[gsg], static_cast(cell_num_nodes)); + } + + for (unsigned int m = 0; m < num_moments_; ++m) + { + const double wn_d2m = d2m_row[m]; + for (std::size_t i = 0; i < cell_num_nodes; ++i) + { + const auto ir = cell_transport_view.MapDOF(i, m, gs_gi); + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) + destination_phi_[ir + gsg] += wn_d2m * b[gsg](i); + } + } + + if (SaveAngularFluxEnabled()) + { + double* cell_psi_data = + &destination_psi_[discretization_.MapDOFLocal(cell, 0, groupset.psi_uk_man_, 0, 0)]; + + for (std::size_t i = 0; i < cell_num_nodes; ++i) + { + const std::size_t imap = + i * groupset_angle_group_stride_ + direction_num * groupset_group_stride_; + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) + cell_psi_data[imap + gsg] = b[gsg](i); + } + } + + for (std::size_t f = 0; f < cell_num_faces; ++f) + { + if (face_orientations[f] != FaceOrientation::OUTGOING) + continue; + + const auto& face = cell.faces[f]; + const bool is_local_face = cell_transport_view.IsFaceLocal(f); + const bool is_boundary_face = not face.has_neighbor; + const bool is_reflecting_boundary_face = + (is_boundary_face and angle_set.GetBoundaries()[face.neighbor_id]->IsReflecting()); + const auto& int_f_shape_i = IntS_shapeI[f]; + + std::vector* psi_nonlocal_outgoing = nullptr; + if (not is_boundary_face and not is_local_face) + psi_nonlocal_outgoing = &workspace_.outgoing_psi_by_face[f]->data; + + const std::size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f); + for (std::size_t fi = 0; fi < num_face_nodes; ++fi) + { + const int i = cell_mapping.MapFaceNode(f, fi); + + if (is_boundary_face) + { + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) + cell_outflow_view.Add( + f, gs_gi + gsg, wt * face_mu_values[f] * b[gsg](i) * int_f_shape_i(i)); + } + + double* psi = nullptr; + if (is_local_face) + psi = fluds.OutgoingPsi(cell, i, as_ss_idx); + else if (not is_boundary_face) + psi = fluds.NLOutgoingPsi(psi_nonlocal_outgoing, fi, as_ss_idx); + else if (is_reflecting_boundary_face) + { + psi = angle_set.PsiReflected(face.neighbor_id, + direction_num, + cell_local_id, + static_cast(f), + static_cast(fi)); + } + + if (psi != nullptr) + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) + psi[gsg] = b[gsg](i); + } + } + + const auto f0 = 1.0 / fac_diamond_difference; + const auto f1 = f0 - 1.0; + for (std::size_t i = 0; i < cell_num_nodes; ++i) + { + const auto ir = discretization_.MapDOFLocal(cell, i, unknown_manager_, polar_level, gs_gi); + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) + psi_sweep_[ir + gsg] = f0 * b[gsg](i) - f1 * psi_sweep_[ir + gsg]; + } + } + + QueueNonlocalOutgoingPsi(workspace_, *async_comm_); +} + +template +void +CBCSweepChunkRZ::Sweep_FixedN(AngleSet& angle_set) +{ + static_assert(NumNodes >= 2 and NumNodes <= 8); + + const auto& groupset = groupset_; + const auto gs_size = groupset.GetNumGroups(); + const auto gs_gi = groupset.first_group; + const auto num_angles_in_as = angle_set.GetNumAngles(); + const auto group_angle_stride = gs_size * num_angles_in_as; + const auto& cell = *cell_; + const auto cell_local_id = cell.local_id; + const auto& cell_mapping = discretization_.GetCellMapping(cell); + const auto& cell_transport_view = cell_transport_views_[cell_local_id]; + auto& cell_outflow_view = cell_outflow_views_[cell_local_id]; + const std::size_t cell_num_faces = cell.faces.size(); + const auto& unit_mats = unit_cell_matrices_[cell_local_id]; + auto& fluds = *fluds_; + const auto& common_data = fluds.GetCommonData(); + const auto& G = unit_mats.intV_shapeI_gradshapeJ; + const auto& M = unit_mats.intV_shapeI_shapeJ; + const auto& M_surf = unit_mats.intS_shapeI_shapeJ; + const auto& IntS_shapeI = unit_mats.intS_shapeI; + const auto& Maux = secondary_unit_cell_matrices_[cell_local_id].intV_shapeI_shapeJ; + const auto& m2d_op = groupset.quadrature->GetMomentToDiscreteOperator(); + const auto& d2m_op = groupset.quadrature->GetDiscreteToMomentOperator(); + const auto& diamond_difference_factor = curvilinear_quadrature_.GetDiamondDifferenceFactor(); + const auto& streaming_operator_factor = curvilinear_quadrature_.GetStreamingOperatorFactor(); + const auto& face_orientations = angle_set.GetSPDS().GetCellFaceOrientations()[cell_local_id]; + const auto& sigma_t = xs_.at(cell.block_id)->GetSigmaTotal(); + + constexpr std::size_t matrix_size = static_cast(NumNodes) * NumNodes; + constexpr auto idx = [](std::size_t i, std::size_t j) -> std::size_t { return i * NumNodes + j; }; + + std::array mass_matrix{}; + std::array aux_matrix{}; + PRAGMA_UNROLL + for (std::size_t i = 0; i < NumNodes; ++i) + { + PRAGMA_UNROLL + for (std::size_t j = 0; j < NumNodes; ++j) + { + mass_matrix[idx(i, j)] = M(i, j); + aux_matrix[idx(i, j)] = Maux(i, j); + } + } + + auto& b = workspace_.rhs_values; + b.resize(static_cast(gs_size) * NumNodes); + auto& sigma_block = workspace_.sigma_t_group_block; + sigma_block.reserve(group_block_size_); + std::array source{}; + auto& moment_dof_map = workspace_.moment_dof_map; + moment_dof_map.resize(static_cast(num_moments_) * NumNodes); + for (unsigned int m = 0; m < num_moments_; ++m) + { + PRAGMA_UNROLL + for (std::size_t i = 0; i < NumNodes; ++i) + moment_dof_map[static_cast(m) * NumNodes + i] = + cell_transport_view.MapDOF(i, m, gs_gi); + } + + auto& face_mu_values = workspace_.face_mu_values; + face_mu_values.assign(cell_num_faces, 0.0); + PrepareNonlocalOutgoingPsi(workspace_, + fluds, + cell, + cell_local_id, + cell_mapping, + cell_transport_view, + group_angle_stride, + face_orientations); + + const auto& as_angle_indices = angle_set.GetAngleIndices(); + for (std::size_t as_ss_idx = 0; as_ss_idx < num_angles_in_as; ++as_ss_idx) + { + const auto direction_num = as_angle_indices[as_ss_idx]; + const auto& omega = groupset.quadrature->GetOmega(direction_num); + const auto wt = groupset.quadrature->GetWeight(direction_num); + + const auto polar_level = direction_polar_level_[direction_num]; + const auto fac_diamond_difference = diamond_difference_factor[direction_num]; + const auto fac_streaming_operator = streaming_operator_factor[direction_num]; + + std::fill(b.begin(), b.end(), 0.0); + + std::array Amat{}; + PRAGMA_UNROLL + for (std::size_t i = 0; i < NumNodes; ++i) + { + PRAGMA_UNROLL + for (std::size_t j = 0; j < NumNodes; ++j) + { + const auto jr = discretization_.MapDOFLocal(cell, j, unknown_manager_, polar_level, gs_gi); + const double streaming = fac_streaming_operator * aux_matrix[idx(i, j)]; + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) + b[gsg * NumNodes + i] += streaming * psi_sweep_[jr + gsg]; + + Amat[idx(i, j)] = omega.Dot(G(i, j)) + streaming; + } + } + + for (std::size_t f = 0; f < cell_num_faces; ++f) + face_mu_values[f] = omega.Dot(cell.faces[f].normal); + + for (std::size_t f = 0; f < cell_num_faces; ++f) + { + if (face_orientations[f] != FaceOrientation::INCOMING) + continue; + + const auto& face = cell.faces[f]; + const bool is_boundary_face = not face.has_neighbor; + const bool is_local_face = not is_boundary_face and cell_transport_view.IsFaceLocal(f); + const auto* face_nodal_mapping = + is_boundary_face + ? nullptr + : &common_data.GetFaceNodalMapping(cell_local_id, static_cast(f)); + const auto incoming_nonlocal_slot = + (is_boundary_face or is_local_face) + ? CBC_FLUDSCommonData::INVALID_FACE_SLOT + : common_data.IncomingFaceSlot(cell_local_id, static_cast(f)); + + const std::size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f); + for (std::size_t fi = 0; fi < num_face_nodes; ++fi) + { + const int i = cell_mapping.MapFaceNode(f, fi); + for (std::size_t fj = 0; fj < num_face_nodes; ++fj) + { + const int j = cell_mapping.MapFaceNode(f, fj); + const double mu_Nij = -face_mu_values[f] * M_surf[f](i, j); + Amat[idx(i, j)] += mu_Nij; + + const double* psi = nullptr; + if (is_local_face) + { + psi = fluds.UpwindPsi(*cell_transport_view.FaceNeighbor(f), + face_nodal_mapping->cell_node_mapping_[fj], + as_ss_idx); + } + else if (not is_boundary_face) + { + psi = fluds.NLUpwindPsi( + incoming_nonlocal_slot, face_nodal_mapping->face_node_mapping_[fj], as_ss_idx); + } + else + { + const bool incident_on_symmetric_boundary = + (face.normal.Dot(normal_vector_boundary_) < -0.999999); + if (not incident_on_symmetric_boundary) + { + psi = angle_set.PsiBoundary(face.neighbor_id, + direction_num, + cell_local_id, + static_cast(f), + static_cast(fj), + gs_gi, + IsSurfaceSourceActive()); + } + } + + if (psi != nullptr) + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) + b[gsg * NumNodes + i] += psi[gsg] * mu_Nij; + } + } + } + + const auto dir_moment_offset = + static_cast(direction_num) * static_cast(num_moments_); + const double* __restrict m2d_row = m2d_op.data() + dir_moment_offset; + const double* __restrict d2m_row = d2m_op.data() + dir_moment_offset; + + for (unsigned int g0 = 0; g0 < gs_size; g0 += group_block_size_) + { + const auto g1 = std::min(g0 + group_block_size_, gs_size); + const auto block_len = g1 - g0; + sigma_block.resize(block_len); + + for (unsigned int gsg = g0; gsg < g1; ++gsg) + { + const auto rel = gsg - g0; + sigma_block[rel] = sigma_t[gs_gi + gsg]; + + source.fill(0.0); + for (unsigned int m = 0; m < num_moments_; ++m) + { + const double w = m2d_row[m]; + const auto moment_offset = static_cast(m) * NumNodes; + PRAGMA_UNROLL + for (std::size_t i = 0; i < NumNodes; ++i) + source[i] += w * source_moments_[moment_dof_map[moment_offset + i] + gsg]; + } + + auto* __restrict bg = &b[static_cast(gsg) * NumNodes]; + PRAGMA_UNROLL + for (std::size_t i = 0; i < NumNodes; ++i) + { + double temp = 0.0; + PRAGMA_UNROLL + for (std::size_t j = 0; j < NumNodes; ++j) + temp += mass_matrix[idx(i, j)] * source[j]; + bg[i] += temp; + } + } + + std::size_t k = 0; +#if __AVX512F__ + for (; k + simd_width <= block_len; k += simd_width) + detail::SimdBatchSolve( + Amat.data(), mass_matrix.data(), &sigma_block[k], &b[(g0 + k) * NumNodes]); +#elif __AVX2__ + for (; k + simd_width <= block_len; k += simd_width) + detail::SimdBatchSolve( + Amat.data(), mass_matrix.data(), &sigma_block[k], &b[(g0 + k) * NumNodes]); +#endif + + for (; k < block_len; ++k) + { + const std::size_t gsg = g0 + k; + const double sigma_tg = sigma_block[k]; + std::array A{}; + PRAGMA_UNROLL + for (std::size_t i = 0; i < NumNodes; ++i) + { + PRAGMA_UNROLL + for (std::size_t j = 0; j < NumNodes; ++j) + A[idx(i, j)] = Amat[idx(i, j)] + sigma_tg * mass_matrix[idx(i, j)]; + } + + auto* __restrict bg = &b[gsg * NumNodes]; + for (std::size_t pivot = 0; pivot < NumNodes; ++pivot) + { + const double inv = 1.0 / A[idx(pivot, pivot)]; + for (std::size_t row = pivot + 1; row < NumNodes; ++row) + { + const double factor = A[idx(row, pivot)] * inv; + bg[row] -= factor * bg[pivot]; + PRAGMA_UNROLL + for (std::size_t col = pivot + 1; col < NumNodes; ++col) + A[idx(row, col)] -= factor * A[idx(pivot, col)]; + } + } + + for (std::size_t pivot = NumNodes; pivot-- > 0;) + { + PRAGMA_UNROLL + for (std::size_t col = pivot + 1; col < NumNodes; ++col) + bg[pivot] -= A[idx(pivot, col)] * bg[col]; + bg[pivot] /= A[idx(pivot, pivot)]; + } + } + + for (unsigned int gsg = g0; gsg < g1; ++gsg) + { + const auto* __restrict bg = &b[static_cast(gsg) * NumNodes]; + for (unsigned int m = 0; m < num_moments_; ++m) + { + const double wn_d2m = d2m_row[m]; + const auto moment_offset = static_cast(m) * NumNodes; + PRAGMA_UNROLL + for (std::size_t i = 0; i < NumNodes; ++i) + destination_phi_[moment_dof_map[moment_offset + i] + gsg] += wn_d2m * bg[i]; + } + } + } + + if (SaveAngularFluxEnabled()) + { + double* cell_psi_data = + &destination_psi_[discretization_.MapDOFLocal(cell, 0, groupset.psi_uk_man_, 0, 0)]; + + PRAGMA_UNROLL + for (std::size_t i = 0; i < NumNodes; ++i) + { + const std::size_t imap = + i * groupset_angle_group_stride_ + direction_num * groupset_group_stride_; + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) + cell_psi_data[imap + gsg] = b[gsg * NumNodes + i]; + } + } + + for (std::size_t f = 0; f < cell_num_faces; ++f) + { + if (face_orientations[f] != FaceOrientation::OUTGOING) + continue; + + const auto& face = cell.faces[f]; + const bool is_local_face = cell_transport_view.IsFaceLocal(f); + const bool is_boundary_face = not face.has_neighbor; + const bool is_reflecting_boundary_face = + (is_boundary_face and angle_set.GetBoundaries()[face.neighbor_id]->IsReflecting()); + const auto& int_f_shape_i = IntS_shapeI[f]; + + std::vector* psi_nonlocal_outgoing = nullptr; + if (not is_boundary_face and not is_local_face) + psi_nonlocal_outgoing = &workspace_.outgoing_psi_by_face[f]->data; + + const std::size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f); + for (std::size_t fi = 0; fi < num_face_nodes; ++fi) + { + const int i = cell_mapping.MapFaceNode(f, fi); + + if (is_boundary_face) + { + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) + cell_outflow_view.Add( + f, gs_gi + gsg, wt * face_mu_values[f] * b[gsg * NumNodes + i] * int_f_shape_i(i)); + } + + double* psi = nullptr; + if (is_local_face) + psi = fluds.OutgoingPsi(cell, i, as_ss_idx); + else if (not is_boundary_face) + psi = fluds.NLOutgoingPsi(psi_nonlocal_outgoing, fi, as_ss_idx); + else if (is_reflecting_boundary_face) + { + psi = angle_set.PsiReflected(face.neighbor_id, + direction_num, + cell_local_id, + static_cast(f), + static_cast(fi)); + } + + if (psi != nullptr) + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) + psi[gsg] = b[gsg * NumNodes + i]; + } + } + + const auto f0 = 1.0 / fac_diamond_difference; + const auto f1 = f0 - 1.0; + PRAGMA_UNROLL + for (std::size_t i = 0; i < NumNodes; ++i) + { + const auto ir = discretization_.MapDOFLocal(cell, i, unknown_manager_, polar_level, gs_gi); + for (std::size_t gsg = 0; gsg < gs_size; ++gsg) + psi_sweep_[ir + gsg] = f0 * b[gsg * NumNodes + i] - f1 * psi_sweep_[ir + gsg]; + } + } + + QueueNonlocalOutgoingPsi(workspace_, *async_comm_); +} + +template void CBCSweepChunkRZ::Sweep_FixedN<2>(AngleSet&); +template void CBCSweepChunkRZ::Sweep_FixedN<3>(AngleSet&); +template void CBCSweepChunkRZ::Sweep_FixedN<4>(AngleSet&); +template void CBCSweepChunkRZ::Sweep_FixedN<5>(AngleSet&); +template void CBCSweepChunkRZ::Sweep_FixedN<6>(AngleSet&); +template void CBCSweepChunkRZ::Sweep_FixedN<7>(AngleSet&); +template void CBCSweepChunkRZ::Sweep_FixedN<8>(AngleSet&); + +} // namespace opensn diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/sweep_chunks/cbc_sweep_chunk_rz.h b/modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/sweep_chunks/cbc_sweep_chunk_rz.h new file mode 100644 index 000000000..6e306299d --- /dev/null +++ b/modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/sweep_chunks/cbc_sweep_chunk_rz.h @@ -0,0 +1,58 @@ +// SPDX-FileCopyrightText: 2026 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#pragma once + +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep_chunks/cbc_sweep_kernels.h" +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep_chunks/sweep_chunk.h" +#include + +namespace opensn +{ + +class DiscreteOrdinatesProblem; +class LBSGroupset; +class CurvilinearProductQuadrature; + +/// Host CBC sweep chunk for 2D RZ curvilinear coordinates. +class CBCSweepChunkRZ : public SweepChunk +{ +public: + CBCSweepChunkRZ(DiscreteOrdinatesProblem& problem, LBSGroupset& groupset); + + void SetAngleSet(AngleSet& angle_set) override; + + void Sweep(AngleSet& angle_set) override; + +private: + using SweepFunc = void (CBCSweepChunkRZ::*)(AngleSet&); + + void Sweep_Generic(AngleSet& angle_set); + + template + void Sweep_FixedN(AngleSet& angle_set); + + /// Curvilinear angular quadrature used for RZ streaming factors. + const CurvilinearProductQuadrature& curvilinear_quadrature_; + /// Secondary spatial discretization cell matrices. + const std::vector& secondary_unit_cell_matrices_; + /// Unknown manager for polar-level sweep dependency storage. + UnknownManager unknown_manager_; + /// Sweeping dependency angular intensity for each polar level. + std::vector psi_sweep_; + /// Direction index to polar level mapping. + std::vector direction_polar_level_; + /// Normal vector used to identify the point/axis of symmetry. + Vector3 normal_vector_boundary_; + + CBC_FLUDS* fluds_ = nullptr; + CBC_AsynchronousCommunicator* async_comm_ = nullptr; + + /// Number of groups solved in one block. + unsigned int group_block_size_; + CBCSweepWorkspace workspace_; + /// Selected sweep implementation. + SweepFunc sweep_impl_ = nullptr; +}; + +} // namespace opensn diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/angle_set/cbc_angle_set.cc b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/angle_set/cbc_angle_set.cc index a90afcc97..e7fef75a8 100644 --- a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/angle_set/cbc_angle_set.cc +++ b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/angle_set/cbc_angle_set.cc @@ -68,7 +68,20 @@ CBC_AngleSet::AngleSetAdvance(SweepChunk& sweep_chunk, AngleSetStatus permission return AngleSetStatus::RECEIVING; if (permission != AngleSetStatus::EXECUTE) - return AngleSetStatus::READY_TO_EXECUTE; + { + const bool all_tasks_completed = (num_completed_tasks_ == task_list_->size()); + const bool all_messages_sent = + (not async_comm_.HasPendingCommunication()) or async_comm_.SendData(); + if (all_tasks_completed and all_messages_sent) + { + for (auto& angleset : following_angle_sets_) + angleset->DecrementCounter(); + executed_ = true; + return AngleSetStatus::FINISHED; + } + + return ready_tasks_.empty() ? AngleSetStatus::NOT_FINISHED : AngleSetStatus::READY_TO_EXECUTE; + } while (not ready_tasks_.empty()) { diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/scheduler/sweep_scheduler.cc b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/scheduler/sweep_scheduler.cc index 5aeaf3aca..be102bf3f 100644 --- a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/scheduler/sweep_scheduler.cc +++ b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/scheduler/sweep_scheduler.cc @@ -11,6 +11,7 @@ #include "framework/runtime.h" #include #include +#include namespace opensn { @@ -18,6 +19,26 @@ namespace opensn namespace { +bool +IsCylindrical(const AngleAggregation& angle_agg) +{ + return angle_agg.GetQuadrature()->GetDimension() == 2 && + angle_agg.GetCoordinateSystem() == CoordinateSystemType::CYLINDRICAL; +} + +bool +HasCurvilinearQuadrature(const AngleAggregation& angle_agg) +{ + return dynamic_cast(angle_agg.GetQuadrature().get()) != + nullptr; +} + +bool +NeedsRZAngleSetDependencies(const AngleAggregation& angle_agg) +{ + return IsCylindrical(angle_agg) && HasCurvilinearQuadrature(angle_agg); +} + bool Compare(const RuleValues& a, const RuleValues& b) { @@ -46,6 +67,60 @@ CompareCylindrical(const RuleValues& a, const RuleValues& b) return a.set_index < b.set_index; } +template +void +ScheduleAngleSets(std::size_t num_angle_sets, + const GetAngleSet& get_angle_set, + SweepChunk& sweep_chunk, + const std::unordered_map>* preceding_angle_sets) +{ + bool finished = false; + std::unordered_set completed; + completed.reserve(num_angle_sets); + while (not finished) + { + finished = true; + for (std::size_t i = 0; i < num_angle_sets; ++i) + { + auto* angle_set = get_angle_set(i); + + if (preceding_angle_sets != nullptr) + { + const auto dep_it = preceding_angle_sets->find(angle_set); + if (dep_it != preceding_angle_sets->end()) + { + bool deps_ready = true; + for (auto* dep : dep_it->second) + { + if (completed.find(dep) == completed.end()) + { + deps_ready = false; + break; + } + } + + if (not deps_ready) + { + const auto status = + angle_set->AngleSetAdvance(sweep_chunk, AngleSetStatus::READY_TO_EXECUTE); + if (status == AngleSetStatus::FINISHED) + completed.insert(angle_set); + else + finished = false; + continue; + } + } + } + + const auto status = angle_set->AngleSetAdvance(sweep_chunk, AngleSetStatus::EXECUTE); + if (status == AngleSetStatus::FINISHED) + completed.insert(angle_set); + if (status != AngleSetStatus::FINISHED) + finished = false; + } + } +} + } // namespace SweepScheduler::SweepScheduler(SchedulingAlgorithm scheduler_type, @@ -53,6 +128,8 @@ SweepScheduler::SweepScheduler(SchedulingAlgorithm scheduler_type, SweepChunk& sweep_chunk) : scheduler_type_(scheduler_type), angle_agg_(angle_agg), sweep_chunk_(sweep_chunk) { + const bool needs_rz_angle_dependencies = NeedsRZAngleSetDependencies(angle_agg_); + if (scheduler_type_ == SchedulingAlgorithm::DEPTH_OF_GRAPH) InitializeAlgoDOG(); @@ -66,18 +143,12 @@ SweepScheduler::SweepScheduler(SchedulingAlgorithm scheduler_type, for (auto& angset : angle_agg_) angset->InitializeDelayedUpstreamData(); - if (scheduler_type_ == SchedulingAlgorithm::DEPTH_OF_GRAPH) + if (needs_rz_angle_dependencies) { - const auto* curvi_quad = - dynamic_cast(angle_agg_.GetQuadrature().get()); - if (curvi_quad && angle_agg_.GetQuadrature()->GetDimension() == 2 && - angle_agg_.GetCoordinateSystem() == CoordinateSystemType::CYLINDRICAL) - { - const auto& following_map = angle_agg_.GetFollowingAngleSetsMap(); - for (const auto& [from, to_set] : following_map) - for (auto* to : to_set) - preceding_angle_sets_[to].insert(from); - } + const auto& following_map = angle_agg_.GetFollowingAngleSetsMap(); + for (const auto& [from, to_set] : following_map) + for (auto* to : to_set) + preceding_angle_sets_[to].insert(from); } // Get local max num messages accross anglesets @@ -97,8 +168,7 @@ SweepScheduler::SweepScheduler(SchedulingAlgorithm scheduler_type, void SweepScheduler::InitializeAlgoDOG() { - const bool is_cylindrical = angle_agg_.GetQuadrature()->GetDimension() == 2 && - angle_agg_.GetCoordinateSystem() == CoordinateSystemType::CYLINDRICAL; + const bool is_cylindrical = IsCylindrical(angle_agg_); std::unordered_map angle_order; const auto* curvi_quad = @@ -113,8 +183,8 @@ SweepScheduler::InitializeAlgoDOG() angle_order.emplace(dir_id, order++); } // Load all anglesets in preperation for sorting - size_t num_anglesets = angle_agg_.GetNumAngleSets(); - for (size_t as = 0; as < num_anglesets; ++as) + std::size_t num_anglesets = angle_agg_.GetNumAngleSets(); + for (std::size_t as = 0; as < num_anglesets; ++as) { auto angleset = angle_agg_[as]; const auto& spds = dynamic_cast(angleset->GetSPDS()); @@ -123,9 +193,9 @@ SweepScheduler::InitializeAlgoDOG() // Find location depth int loc_depth = -1; - for (size_t level = 0; level < leveled_graph.size(); ++level) + for (std::size_t level = 0; level < leveled_graph.size(); ++level) { - for (size_t index = 0; index < leveled_graph[level].item_id.size(); ++index) + for (std::size_t index = 0; index < leveled_graph[level].item_id.size(); ++index) { if (leveled_graph[level].item_id[index] == opensn::mpi_comm.rank()) { @@ -171,18 +241,13 @@ SweepScheduler::ScheduleAlgoDOG(SweepChunk& sweep_chunk) for (auto& angle_set : angle_agg_) angle_set->ResetDependencyCounter(); - bool finished = false; - while (not finished) - { - finished = true; - for (auto& rule_value : rule_values_) - { - auto angleset = rule_value.angle_set; - AngleSetStatus status = angleset->AngleSetAdvance(sweep_chunk, AngleSetStatus::EXECUTE); - if (status != AngleSetStatus::FINISHED) - finished = false; - } - } + const auto* const preceding_angle_sets = + NeedsRZAngleSetDependencies(angle_agg_) ? &preceding_angle_sets_ : nullptr; + ScheduleAngleSets( + rule_values_.size(), + [this](const std::size_t i) { return rule_values_[i].angle_set.get(); }, + sweep_chunk, + preceding_angle_sets); // Receive delayed data opensn::mpi_comm.barrier(); @@ -213,19 +278,13 @@ SweepScheduler::ScheduleAlgoFIFO(SweepChunk& sweep_chunk) for (auto& angle_set : angle_agg_) angle_set->ResetDependencyCounter(); - // Loop over AngleSetGroups - bool finished = false; - while (not finished) - { - finished = true; - - for (auto& angle_set : angle_agg_) - { - AngleSetStatus status = angle_set->AngleSetAdvance(sweep_chunk, AngleSetStatus::EXECUTE); - if (status != AngleSetStatus::FINISHED) - finished = false; - } // for angleset - } // while not finished + const auto* const preceding_angle_sets = + NeedsRZAngleSetDependencies(angle_agg_) ? &preceding_angle_sets_ : nullptr; + ScheduleAngleSets( + angle_agg_.GetNumAngleSets(), + [this](const std::size_t i) { return angle_agg_[i].get(); }, + sweep_chunk, + preceding_angle_sets); // Receive delayed data opensn::mpi_comm.barrier(); diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/scheduler/sweep_scheduler.h b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/scheduler/sweep_scheduler.h index 9cae0b817..46ad85877 100644 --- a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/scheduler/sweep_scheduler.h +++ b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/scheduler/sweep_scheduler.h @@ -6,6 +6,7 @@ #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/angle_aggregation/angle_aggregation.h" #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep_chunks/sweep_chunk.h" #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/scheduler/spmd_threadpool.h" +#include #include #include @@ -30,7 +31,7 @@ struct RuleValues int sign_of_omegay; int sign_of_omegaz; int azimuthal_order; - size_t set_index; + std::size_t set_index; explicit RuleValues(std::shared_ptr& ref_as) : angle_set(ref_as), @@ -78,7 +79,7 @@ class SweepScheduler std::vector rule_values_; - /// Angle set dependencies (preceding sets) used by DOG scheduling. + /// Angle set dependencies (preceding sets) used by RZ scheduling. std::unordered_map> preceding_angle_sets_; SPMD_ThreadPool pool_; diff --git a/test/python/modules/linear_boltzmann_solvers/transport_keigen_cyl/tests.json b/test/python/modules/linear_boltzmann_solvers/transport_keigen_cyl/tests.json index 84681d4ff..f2c551f7c 100644 --- a/test/python/modules/linear_boltzmann_solvers/transport_keigen_cyl/tests.json +++ b/test/python/modules/linear_boltzmann_solvers/transport_keigen_cyl/tests.json @@ -12,6 +12,19 @@ } ] }, + { + "file": "transport_2d_cyl_disk_keigen_vacuum_bcs_cbc.py", + "comment": "RZ disk k-eigenvalue vacuum boundaries (3D disk reference) (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "KeyValuePair", + "key": "K_EIGEN", + "goldvalue": 0.0900968, + "abs_tol": 0.005 + } + ] + }, { "file": "transport_2d_cyl_disk_keigen_vacuum_gmsh_struct.py", "comment": "RZ disk k-eigenvalue vacuum boundaries (gmsh structured)", @@ -25,6 +38,19 @@ } ] }, + { + "file": "transport_2d_cyl_disk_keigen_vacuum_gmsh_struct_cbc.py", + "comment": "RZ disk k-eigenvalue vacuum boundaries (gmsh structured) (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "KeyValuePair", + "key": "K_EIGEN", + "goldvalue": 0.0900968, + "abs_tol": 0.005 + } + ] + }, { "file": "transport_2d_cyl_disk_keigen_vacuum_gmsh_unstruct.py", "comment": "RZ disk k-eigenvalue vacuum boundaries (gmsh unstructured)", @@ -37,5 +63,18 @@ "abs_tol": 0.005 } ] + }, + { + "file": "transport_2d_cyl_disk_keigen_vacuum_gmsh_unstruct_cbc.py", + "comment": "RZ disk k-eigenvalue vacuum boundaries (gmsh unstructured) (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "KeyValuePair", + "key": "K_EIGEN", + "goldvalue": 0.0900968, + "abs_tol": 0.005 + } + ] } ] diff --git a/test/python/modules/linear_boltzmann_solvers/transport_keigen_cyl/transport_2d_cyl_disk_keigen_vacuum_bcs_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_keigen_cyl/transport_2d_cyl_disk_keigen_vacuum_bcs_cbc.py new file mode 100644 index 000000000..5139dc4f7 --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_keigen_cyl/transport_2d_cyl_disk_keigen_vacuum_bcs_cbc.py @@ -0,0 +1,101 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +RZ k-eigenvalue test with vacuum boundaries. + +Cylinder with r=0.2 m and h=0.1 m, axisymmetric (rmin reflecting), vacuum +boundaries at rmax/zmin/zmax. +K_EIGEN = 0.0900968 from 3D disk keigen vacuum transport solution +(transport_3d_disk_keigen_vacuum_bcs.py). +""" + +import os +import sys +import math + +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.mesh import OrthogonalMeshGenerator + from pyopensn.xs import MultiGroupXS + from pyopensn.aquad import GLCProductQuadrature2DRZ + from pyopensn.solver import DiscreteOrdinatesCurvilinearProblem, PowerIterationKEigenSolver + + +if __name__ == "__main__": + rmin = 0.0 + rmax = 0.2 + zmin = 0.0 + zmax = 0.1 + nr = 100 + nz = 50 + r_nodes = [rmin + i * ((rmax - rmin) / nr) for i in range(nr + 1)] + z_nodes = [zmin + i * ((zmax - zmin) / nz) for i in range(nz + 1)] + meshgen = OrthogonalMeshGenerator(node_sets=[r_nodes, z_nodes], coord_sys="cylindrical") + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + + xs = MultiGroupXS() + xs.LoadFromOpenSn("../../../../assets/xs/simple_fissile_1g.xs") + + quad = GLCProductQuadrature2DRZ(n_polar=8, n_azimuthal=8, scattering_order=0) + bcs = [ + {"name": "rmin", "type": "reflecting"}, # axis symmetry + {"name": "rmax", "type": "vacuum"}, + {"name": "zmin", "type": "vacuum"}, + {"name": "zmax", "type": "vacuum"}, + ] + + phys = DiscreteOrdinatesCurvilinearProblem( + mesh=grid, + sweep_type="CBC", + num_groups=1, + groupsets=[ + { + "groups_from_to": (0, 0), + "angular_quadrature": quad, + "angle_aggregation_type": "single", + "inner_linear_method": "petsc_gmres", + "l_max_its": 200, + "l_abs_tol": 1.0e-10, + "gmres_restart_interval": 100, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + boundary_conditions=bcs, + options={ + "use_precursors": False, + "verbose_inner_iterations": False, + }, + ) + + k_solver = PowerIterationKEigenSolver(problem=phys, k_tol=1.0e-8) + k_solver.Initialize() + k_solver.Execute() + + k = k_solver.GetEigenvalue() + + sigma_t = xs.sigma_t[0] + sigma_a = xs.sigma_a[0] if len(xs.sigma_a) > 0 else sigma_t + nu_sigma_f = xs.nu_sigma_f[0] + + D = 1.0 / (3.0 * sigma_t) + R = rmax - rmin + H = zmax - zmin + z_ex = 0.7104 * D + R_e = R + z_ex + H_e = H + 2.0 * z_ex + B2 = (2.405 / R_e) ** 2 + (math.pi / H_e) ** 2 + k_ref = nu_sigma_f / (sigma_a + D * B2) + + tol = float(os.environ.get("OPENSN_TOL", "5e-2")) + ok = abs(k - k_ref) <= tol + + if rank == 0: + print(f"K_EIGEN {k:.12e}") + print(f"K_REF {k_ref:.12e}") + print(f"TOL {tol:.6e}") + print(f"PASS {int(ok)}") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_keigen_cyl/transport_2d_cyl_disk_keigen_vacuum_gmsh_struct_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_keigen_cyl/transport_2d_cyl_disk_keigen_vacuum_gmsh_struct_cbc.py new file mode 100644 index 000000000..6e7087095 --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_keigen_cyl/transport_2d_cyl_disk_keigen_vacuum_gmsh_struct_cbc.py @@ -0,0 +1,99 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +RZ k-eigenvalue test with vacuum boundaries and gmsh orthogonal mesh. + +Cylinder with r=0.2 m and h=0.1 m, axisymmetric (rmin reflecting), vacuum +boundaries at rmax/zmin/zmax. +K_EIGEN = 0.0900968 from 3D disk keigen vacuum transport solution +(transport_3d_disk_keigen_vacuum_bcs.py). +""" + +import os +import sys +import math + +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.mesh import FromFileMeshGenerator + from pyopensn.xs import MultiGroupXS + from pyopensn.aquad import GLCProductQuadrature2DRZ + from pyopensn.solver import DiscreteOrdinatesCurvilinearProblem, PowerIterationKEigenSolver + + +if __name__ == "__main__": + mesh_file = "../../../../assets/mesh/rz_disk_struct_quad.msh" + meshgen = FromFileMeshGenerator(filename=mesh_file, coord_sys="cylindrical") + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + + rmin = 0.0 + rmax = 0.2 + zmin = 0.0 + zmax = 0.1 + + xs = MultiGroupXS() + xs.LoadFromOpenSn("../../../../assets/xs/simple_fissile_1g.xs") + + quad = GLCProductQuadrature2DRZ(n_polar=16, n_azimuthal=32, scattering_order=0) + bcs = [ + {"name": "rmin", "type": "reflecting"}, + {"name": "rmax", "type": "vacuum"}, + {"name": "zmin", "type": "vacuum"}, + {"name": "zmax", "type": "vacuum"}, + ] + + phys = DiscreteOrdinatesCurvilinearProblem( + mesh=grid, + sweep_type="CBC", + num_groups=1, + groupsets=[ + { + "groups_from_to": (0, 0), + "angular_quadrature": quad, + "angle_aggregation_type": "single", + "inner_linear_method": "petsc_gmres", + "l_max_its": 200, + "l_abs_tol": 1.0e-10, + "gmres_restart_interval": 100, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + boundary_conditions=bcs, + options={ + "use_precursors": False, + "verbose_inner_iterations": False, + }, + ) + + k_solver = PowerIterationKEigenSolver(problem=phys, k_tol=1.0e-8) + k_solver.Initialize() + k_solver.Execute() + + k = k_solver.GetEigenvalue() + + sigma_t = xs.sigma_t[0] + sigma_a = xs.sigma_a[0] if len(xs.sigma_a) > 0 else sigma_t + nu_sigma_f = xs.nu_sigma_f[0] + + D = 1.0 / (3.0 * sigma_t) + R = rmax - rmin + H = zmax - zmin + z_ex = 0.7104 * D + R_e = R + z_ex + H_e = H + 2.0 * z_ex + B2 = (2.405 / R_e) ** 2 + (math.pi / H_e) ** 2 + k_ref = nu_sigma_f / (sigma_a + D * B2) + + tol = float(os.environ.get("OPENSN_TOL", "5e-2")) + ok = abs(k - k_ref) <= tol + + if rank == 0: + print(f"K_EIGEN {k:.12e}") + print(f"K_REF {k_ref:.12e}") + print(f"TOL {tol:.6e}") + print(f"PASS {int(ok)}") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_keigen_cyl/transport_2d_cyl_disk_keigen_vacuum_gmsh_unstruct_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_keigen_cyl/transport_2d_cyl_disk_keigen_vacuum_gmsh_unstruct_cbc.py new file mode 100644 index 000000000..1b0a30046 --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_keigen_cyl/transport_2d_cyl_disk_keigen_vacuum_gmsh_unstruct_cbc.py @@ -0,0 +1,99 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +RZ k-eigenvalue test with vacuum boundaries and gmsh unstructured mesh. + +Cylinder with r=0.2 m and h=0.1 m, axisymmetric (rmin reflecting), vacuum +boundaries at rmax/zmin/zmax. +K_EIGEN = 0.0900968 from 3D disk keigen vacuum transport solution +(transport_3d_disk_keigen_vacuum_bcs.py). +""" + +import os +import sys +import math + +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.mesh import FromFileMeshGenerator + from pyopensn.xs import MultiGroupXS + from pyopensn.aquad import GLCProductQuadrature2DRZ + from pyopensn.solver import DiscreteOrdinatesCurvilinearProblem, PowerIterationKEigenSolver + + +if __name__ == "__main__": + mesh_file = "../../../../assets/mesh/rz_disk_unstructured.msh" + meshgen = FromFileMeshGenerator(filename=mesh_file, coord_sys="cylindrical") + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + + rmin = 0.0 + rmax = 0.2 + zmin = 0.0 + zmax = 0.1 + + xs = MultiGroupXS() + xs.LoadFromOpenSn("../../../../assets/xs/simple_fissile_1g.xs") + + quad = GLCProductQuadrature2DRZ(n_polar=16, n_azimuthal=32, scattering_order=0) + bcs = [ + {"name": "rmin", "type": "reflecting"}, + {"name": "rmax", "type": "vacuum"}, + {"name": "zmin", "type": "vacuum"}, + {"name": "zmax", "type": "vacuum"}, + ] + + phys = DiscreteOrdinatesCurvilinearProblem( + mesh=grid, + sweep_type="CBC", + num_groups=1, + groupsets=[ + { + "groups_from_to": (0, 0), + "angular_quadrature": quad, + "angle_aggregation_type": "single", + "inner_linear_method": "petsc_gmres", + "l_max_its": 200, + "l_abs_tol": 1.0e-10, + "gmres_restart_interval": 100, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + boundary_conditions=bcs, + options={ + "use_precursors": False, + "verbose_inner_iterations": False, + }, + ) + + k_solver = PowerIterationKEigenSolver(problem=phys, k_tol=1.0e-8) + k_solver.Initialize() + k_solver.Execute() + + k = k_solver.GetEigenvalue() + + sigma_t = xs.sigma_t[0] + sigma_a = xs.sigma_a[0] if len(xs.sigma_a) > 0 else sigma_t + nu_sigma_f = xs.nu_sigma_f[0] + + D = 1.0 / (3.0 * sigma_t) + R = rmax - rmin + H = zmax - zmin + z_ex = 0.7104 * D + R_e = R + z_ex + H_e = H + 2.0 * z_ex + B2 = (2.405 / R_e) ** 2 + (math.pi / H_e) ** 2 + k_ref = nu_sigma_f / (sigma_a + D * B2) + + tol = float(os.environ.get("OPENSN_TOL", "5e-2")) + ok = abs(k - k_ref) <= tol + + if rank == 0: + print(f"K_EIGEN {k:.12e}") + print(f"K_REF {k_ref:.12e}") + print(f"TOL {tol:.6e}") + print(f"PASS {int(ok)}") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/tests.json b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/tests.json index f69283064..a2ba021cc 100644 --- a/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/tests.json +++ b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/tests.json @@ -12,6 +12,19 @@ } ] }, + { + "file": "transport_2d_cyl_1g_cbc.py", + "comment": "2D LinearBSolver Cylindrical Test mono-energetic - PWLD (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "KeyValuePair", + "key": "Max-value=", + "goldvalue": 1.0, + "abs_tol": 1e-09 + } + ] + }, { "file": "transport_2d_cyl_2g.py", "comment": "2D LinearBSolver Cylindrical Test multi-group - PWLD", @@ -31,6 +44,25 @@ } ] }, + { + "file": "transport_2d_cyl_2g_cbc.py", + "comment": "2D LinearBSolver Cylindrical Test multi-group - PWLD (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "KeyValuePair", + "key": "Max-valueG1=", + "goldvalue": 1.0, + "abs_tol": 1e-09 + }, + { + "type": "KeyValuePair", + "key": "Max-valueG2=", + "goldvalue": 0.25, + "abs_tol": 1e-09 + } + ] + }, { "file": "transport_2d_cyl_reflective_uniform_1g.py", "comment": "RZ gmsh case 01: reflective 1g uniform", @@ -44,6 +76,19 @@ } ] }, + { + "file": "transport_2d_cyl_reflective_uniform_1g_cbc.py", + "comment": "RZ gmsh case 01: reflective 1g uniform (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "KeyValuePair", + "key": "MAX", + "goldvalue": 0.8437473489223, + "abs_tol": 1e-09 + } + ] + }, { "file": "transport_2d_cyl_reflective_uniform_2g.py", "comment": "RZ gmsh case 02: reflective 2g uniform", @@ -63,6 +108,25 @@ } ] }, + { + "file": "transport_2d_cyl_reflective_uniform_2g_cbc.py", + "comment": "RZ gmsh case 02: reflective 2g uniform (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "KeyValuePair", + "key": "G0_MAX", + "goldvalue": 0.8437473489224, + "abs_tol": 1e-09 + }, + { + "type": "KeyValuePair", + "key": "G1_MAX", + "goldvalue": 0.3658623519312, + "abs_tol": 1e-09 + } + ] + }, { "file": "transport_2d_cyl_annulus_zero_solution.py", "comment": "RZ gmsh case 04: annulus zero-solution", @@ -76,6 +140,19 @@ } ] }, + { + "file": "transport_2d_cyl_annulus_zero_solution_cbc.py", + "comment": "RZ gmsh case 04: annulus zero-solution (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "KeyValuePair", + "key": "MAX", + "goldvalue": 0.0, + "abs_tol": 1e-12 + } + ] + }, { "file": "transport_2d_cyl_unstructured_reflective_uniform_1g.py", "comment": "RZ gmsh case 05: unstructured reflective", @@ -89,6 +166,19 @@ } ] }, + { + "file": "transport_2d_cyl_unstructured_reflective_uniform_1g_cbc.py", + "comment": "RZ gmsh case 05: unstructured reflective (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "KeyValuePair", + "key": "MAX", + "goldvalue": 0.9065113759285, + "abs_tol": 1e-09 + } + ] + }, { "file": "transport_2d_cyl_2g_scatter_vacuum.py", "comment": "RZ gmsh case 07: 2g scatter with vacuum boundaries", @@ -108,6 +198,25 @@ } ] }, + { + "file": "transport_2d_cyl_2g_scatter_vacuum_cbc.py", + "comment": "RZ gmsh case 07: 2g scatter with vacuum boundaries (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "KeyValuePair", + "key": "G0_MAX", + "goldvalue": 0.8429208738494, + "abs_tol": 2e-09 + }, + { + "type": "KeyValuePair", + "key": "G1_MAX", + "goldvalue": 0.3655486193072, + "abs_tol": 3e-09 + } + ] + }, { "file": "transport_2d_cyl_zero_source_vacuum_1g_struct.py", "comment": "RZ analytic: zero source with vacuum boundaries (structured)", @@ -121,6 +230,19 @@ } ] }, + { + "file": "transport_2d_cyl_zero_source_vacuum_1g_struct_cbc.py", + "comment": "RZ analytic: zero source with vacuum boundaries (structured) (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "KeyValuePair", + "key": "PHI_MAX", + "goldvalue": 0.0, + "abs_tol": 1e-12 + } + ] + }, { "file": "transport_2d_cyl_disk_leakage_vacuum_bcs.py", "comment": "RZ disk leakage: vacuum bcs compared to 3D reference", @@ -158,6 +280,43 @@ } ] }, + { + "file": "transport_2d_cyl_disk_leakage_vacuum_bcs_cbc.py", + "comment": "RZ disk leakage: vacuum bcs compared to 3D reference (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "KeyValuePair", + "key": "Top leakage=", + "goldvalue": 36.173790857624596, + "abs_tol": 1e-06 + }, + { + "type": "KeyValuePair", + "key": "Bottom leakage=", + "goldvalue": 36.17379086382096, + "abs_tol": 1e-06 + }, + { + "type": "KeyValuePair", + "key": "Side leakage=", + "goldvalue": 39.74878117223803, + "abs_tol": 1e-06 + }, + { + "type": "KeyValuePair", + "key": "Total leakage=", + "goldvalue": 112.0963628936836, + "abs_tol": 1e-06 + }, + { + "type": "KeyValuePair", + "key": "Leakage/Source=", + "goldvalue": 0.9144851226179025, + "abs_tol": 1e-06 + } + ] + }, { "file": "transport_2d_cyl_disk_leakage_reflecting_bcs.py", "comment": "RZ disk leakage: reflecting z, vacuum side compared to 3D reference", @@ -195,6 +354,43 @@ } ] }, + { + "file": "transport_2d_cyl_disk_leakage_reflecting_bcs_cbc.py", + "comment": "RZ disk leakage: reflecting z, vacuum side compared to 3D reference (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "KeyValuePair", + "key": "Top leakage=", + "goldvalue": -6.190966583965322e-15, + "abs_tol": 1e-12 + }, + { + "type": "KeyValuePair", + "key": "Bottom leakage=", + "goldvalue": -5.231802747012948e-16, + "abs_tol": 1e-12 + }, + { + "type": "KeyValuePair", + "key": "Side leakage=", + "goldvalue": 97.24863271745009, + "abs_tol": 8e-06 + }, + { + "type": "KeyValuePair", + "key": "Total leakage=", + "goldvalue": 97.24863271745009, + "abs_tol": 8e-06 + }, + { + "type": "KeyValuePair", + "key": "Leakage/Source=", + "goldvalue": 0.7933569432523587, + "abs_tol": 1e-06 + } + ] + }, { "file": "transport_2d_cyl_disk_source_vacuum_hdpe_struct.py", "comment": "RZ disk fixed-source vacuum boundaries (HDPE OpenMC XS, PHI_MAX reference)", @@ -208,6 +404,19 @@ } ] }, + { + "file": "transport_2d_cyl_disk_source_vacuum_hdpe_struct_cbc.py", + "comment": "RZ disk fixed-source vacuum boundaries (HDPE OpenMC XS, PHI_MAX reference) (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "KeyValuePair", + "key": "PHI_MAX", + "goldvalue": 0.1224546944901, + "abs_tol": 0.0001 + } + ] + }, { "file": "transport_2d_cyl_2g_scatter_vacuum_ucoarse.py", "comment": "RZ gmsh unstructured 2g scatter vacuum (coarse)", @@ -221,6 +430,19 @@ } ] }, + { + "file": "transport_2d_cyl_2g_scatter_vacuum_ucoarse_cbc.py", + "comment": "RZ gmsh unstructured 2g scatter vacuum (coarse) (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "KeyValuePair", + "key": "PASS", + "goldvalue": 1, + "abs_tol": 1e-12 + } + ] + }, { "file": "transport_2d_cyl_2g_scatter_vacuum_umed.py", "comment": "RZ gmsh unstructured 2g scatter vacuum (medium)", @@ -234,6 +456,19 @@ } ] }, + { + "file": "transport_2d_cyl_2g_scatter_vacuum_umed_cbc.py", + "comment": "RZ gmsh unstructured 2g scatter vacuum (medium) (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "KeyValuePair", + "key": "PASS", + "goldvalue": 1, + "abs_tol": 1e-12 + } + ] + }, { "file": "transport_2d_cyl_2g_scatter_vacuum_ufine.py", "comment": "RZ gmsh unstructured 2g scatter vacuum (fine)", @@ -246,5 +481,18 @@ "abs_tol": 1e-12 } ] + }, + { + "file": "transport_2d_cyl_2g_scatter_vacuum_ufine_cbc.py", + "comment": "RZ gmsh unstructured 2g scatter vacuum (fine) (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "KeyValuePair", + "key": "PASS", + "goldvalue": 1, + "abs_tol": 1e-12 + } + ] } ] diff --git a/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_1g_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_1g_cbc.py new file mode 100644 index 000000000..0cf7670d2 --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_1g_cbc.py @@ -0,0 +1,105 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +transport_2d_cyl_1_monoenergetic.py: 1-group fixed-source with leakage +Expected: Max-value=1.0 +""" + +import os +import sys +import math + +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.mesh import OrthogonalMeshGenerator + from pyopensn.logvol import RPPLogicalVolume + from pyopensn.xs import MultiGroupXS + from pyopensn.source import VolumetricSource + from pyopensn.aquad import GLCProductQuadrature2DRZ + from pyopensn.solver import DiscreteOrdinatesCurvilinearProblem, SteadyStateSourceSolver + from pyopensn.fieldfunc import FieldFunctionInterpolationVolume + +if __name__ == "__main__": + + num_procs = 4 + if size != num_procs: + sys.exit(f"Incorrect number of processors. Expected {num_procs} but got {size}.") + + # Setup mesh + dim = 2 + length = [1.0, 2.0] + ncells = [50, 100] + nodes = [] + for d in range(dim): + delta = length[d] / ncells[d] + node_set = [i * delta for i in range(ncells[d] + 1)] + nodes.append(node_set) + meshgen = OrthogonalMeshGenerator(node_sets=[nodes[0], nodes[1]], coord_sys="cylindrical") + grid = meshgen.Execute() + + # Set block IDs + vol0 = RPPLogicalVolume( + xmin=0.0, + xmax=length[0], + ymin=0.0, + ymax=length[1], + infz=True, + ) + grid.SetBlockIDFromLogicalVolume(vol0, 0, True) + + # Define materials and sources (one-group problem) + ngrp = 1 + sigmat = 25.0 + ratioc = 0.1 + xs_1g = MultiGroupXS() + xs_1g.CreateSimpleOneGroup(sigmat, ratioc) + source = [] + source.append(sigmat * (1 - ratioc)) + mg_src = VolumetricSource(block_ids=[0], group_strength=source) + + # Angular quadrature + pquad = GLCProductQuadrature2DRZ(n_polar=4, n_azimuthal=8, scattering_order=0) + + # Create the curvilinear solver (for cylindrical geometry) + phys = DiscreteOrdinatesCurvilinearProblem( + mesh=grid, + sweep_type="CBC", + num_groups=ngrp, + groupsets=[ + { + "groups_from_to": (0, ngrp - 1), + "angular_quadrature": pquad, + "angle_aggregation_type": "azimuthal", + "inner_linear_method": "petsc_gmres", + "l_max_its": 100, + "l_abs_tol": 1.0e-12, + }, + ], + xs_map=[ + {"block_ids": [0], "xs": xs_1g}, + ], + volumetric_sources=[mg_src], + boundary_conditions=[{"name": "rmin", "type": "reflecting"}], + ) + ss_solver = SteadyStateSourceSolver(problem=phys) + ss_solver.Initialize() + ss_solver.Execute() + + # Field functions + fflist = phys.GetScalarFluxFieldFunction(only_scalar_flux=False) + + # Perform a volume integration over vol0 (compute the maximum value) + ffi1 = FieldFunctionInterpolationVolume() + curffi = ffi1 + curffi.SetOperationType("max") + curffi.SetLogicalVolume(vol0) + curffi.AddFieldFunction(fflist[0][0]) + curffi.Execute() + maxval = curffi.GetValue() + if rank == 0: + print(f"Max-value={maxval:.5f}") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_2g_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_2g_cbc.py new file mode 100644 index 000000000..222e11355 --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_2g_cbc.py @@ -0,0 +1,114 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +transport_2d_cyl_2_multigroup.py: 2-group fixed-source with leakage +Expected: Max-valueG1=1.0, Max-valueG2=0.25. +""" + +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.mesh import OrthogonalMeshGenerator + from pyopensn.logvol import RPPLogicalVolume + from pyopensn.xs import MultiGroupXS + from pyopensn.source import VolumetricSource + from pyopensn.aquad import GLCProductQuadrature2DRZ + from pyopensn.solver import DiscreteOrdinatesCurvilinearProblem, SteadyStateSourceSolver + from pyopensn.fieldfunc import FieldFunctionInterpolationVolume + +if __name__ == "__main__": + + num_procs = 4 + if size != num_procs: + sys.exit(f"Incorrect number of processors. Expected {num_procs} but got {size}.") + + # Setup mesh + dim = 2 + length = [1.0, 2.0] + ncells = [50, 100] + nodes = [] + for d in range(dim): + delta = length[d] / ncells[d] + node_set = [i * delta for i in range(ncells[d] + 1)] + nodes.append(node_set) + meshgen = OrthogonalMeshGenerator(node_sets=[nodes[0], nodes[1]], coord_sys="cylindrical") + grid = meshgen.Execute() + + # Define a logical volume covering the entire domain and set block IDs + vol0 = RPPLogicalVolume(xmin=0.0, xmax=length[0], ymin=0.0, ymax=length[1], infz=True) + grid.SetBlockIDFromLogicalVolume(vol0, 0, True) + + # Define material properties and the volumetric source for a two-group problem + ngrp = 2 + sigmat = 20.0 + ratioc = 0.4 + + # Create source strengths (group 1 nonzero; group 2 zero) + source = [0.0] * ngrp + source[0] = sigmat * (1 - 0.5 * ratioc) + mg_src = VolumetricSource(block_ids=[0], group_strength=source) + + # Load cross-section data from file + xs_data = MultiGroupXS() + xs_data.LoadFromOpenSn("../../../../assets/xs/transport_2d_cyl_2g.xs") + + # Angular quadrature + pquad = GLCProductQuadrature2DRZ(n_polar=4, n_azimuthal=8, scattering_order=0) + + # Create and configure the curvilinear solver for cylindrical geometry + phys = DiscreteOrdinatesCurvilinearProblem( + mesh=grid, + sweep_type="CBC", + num_groups=ngrp, + groupsets=[ + { + "groups_from_to": (0, ngrp - 1), + "angular_quadrature": pquad, + "angle_aggregation_type": "azimuthal", + "inner_linear_method": "petsc_gmres", + "l_max_its": 100, + "l_abs_tol": 1.0e-12, + "apply_wgdsa": True, + "wgdsa_l_abs_tol": 1.0e-9, + "wgdsa_l_max_its": 50, + }, + ], + xs_map=[ + {"block_ids": [0], "xs": xs_data}, + ], + volumetric_sources=[mg_src], + boundary_conditions=[{"name": "rmin", "type": "reflecting"}], + ) + ss_solver = SteadyStateSourceSolver(problem=phys) + ss_solver.Initialize() + ss_solver.Execute() + + # Field functions + fflist = phys.GetScalarFluxFieldFunction(only_scalar_flux=False) + + # Volume integration for energy group 1 + ffi1 = FieldFunctionInterpolationVolume() + ffi1.SetOperationType("max") + ffi1.SetLogicalVolume(vol0) + ffi1.AddFieldFunction(fflist[0][0]) + ffi1.Execute() + maxval = ffi1.GetValue() + if rank == 0: + print(f"Max-valueG1={maxval:.5f}") + + # Volume integration for energy group 2 + ffi2 = FieldFunctionInterpolationVolume() + ffi2.SetOperationType("max") + ffi2.SetLogicalVolume(vol0) + ffi2.AddFieldFunction(fflist[1][0]) + ffi2.Execute() + maxval = ffi2.GetValue() + if rank == 0: + print(f"Max-valueG2={maxval:.5f}") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_2g_scatter_vacuum_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_2g_scatter_vacuum_cbc.py new file mode 100644 index 000000000..43181a07d --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_2g_scatter_vacuum_cbc.py @@ -0,0 +1,91 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +transport_2d_cyl_2g_scatter_vacuum.py: 2-group scatter with vacuum boundaries +Expected: G0_MAX=0.8429208738494, G1_MAX=0.3655486193072 from 3D Cartesian calculation +""" + +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.mesh import OrthogonalMeshGenerator + from pyopensn.logvol import RPPLogicalVolume + from pyopensn.xs import MultiGroupXS + from pyopensn.source import VolumetricSource + from pyopensn.aquad import GLCProductQuadrature2DRZ + from pyopensn.solver import DiscreteOrdinatesCurvilinearProblem, SteadyStateSourceSolver + from pyopensn.fieldfunc import FieldFunctionInterpolationVolume + + +if __name__ == "__main__": + # Structured cylindrical mesh on 0<=r<=1, 0<=z<=2. + nr = 50 + nz = 100 + r_nodes = [i / nr for i in range(nr + 1)] + z_nodes = [2.0 * i / nz for i in range(nz + 1)] + meshgen = OrthogonalMeshGenerator(node_sets=[r_nodes, z_nodes], coord_sys="cylindrical") + grid = meshgen.Execute() + + vol = RPPLogicalVolume(xmin=0.0, xmax=1.0, ymin=0.0, ymax=2.0, infz=True) + grid.SetBlockIDFromLogicalVolume(vol, 0, True) + + xs = MultiGroupXS() + xs.LoadFromOpenSn("../../../../assets/xs/transport_2d_cyl_2g_2.xs") + + src = VolumetricSource(block_ids=[0], group_strength=[1.5, 0.8]) + + quad = GLCProductQuadrature2DRZ(n_polar=4, n_azimuthal=8, scattering_order=0) + problem = DiscreteOrdinatesCurvilinearProblem( + mesh=grid, + sweep_type="CBC", + num_groups=2, + groupsets=[ + { + "groups_from_to": (0, 1), + "angular_quadrature": quad, + "angle_aggregation_type": "azimuthal", + "inner_linear_method": "classic_richardson", + "l_abs_tol": 1.0e-7, + "l_max_its": 5000, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + volumetric_sources=[src], + boundary_conditions=[ + {"name": "rmin", "type": "reflecting"}, + {"name": "rmax", "type": "vacuum"}, + {"name": "zmin", "type": "vacuum"}, + {"name": "zmax", "type": "vacuum"}, + ], + ) + + solver = SteadyStateSourceSolver(problem=problem) + solver.Initialize() + solver.Execute() + + fflist = problem.GetScalarFluxFieldFunction(only_scalar_flux=False) + + ffi_g0 = FieldFunctionInterpolationVolume() + ffi_g0.SetOperationType("max") + ffi_g0.SetLogicalVolume(vol) + ffi_g0.AddFieldFunction(fflist[0][0]) + ffi_g0.Execute() + phi_g0 = ffi_g0.GetValue() + + ffi_g1 = FieldFunctionInterpolationVolume() + ffi_g1.SetOperationType("max") + ffi_g1.SetLogicalVolume(vol) + ffi_g1.AddFieldFunction(fflist[1][0]) + ffi_g1.Execute() + phi_g1 = ffi_g1.GetValue() + + if rank == 0: + print(f"G0_MAX {phi_g0:.12e}") + print(f"G1_MAX {phi_g1:.12e}") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_2g_scatter_vacuum_ucoarse_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_2g_scatter_vacuum_ucoarse_cbc.py new file mode 100644 index 000000000..997e91afc --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_2g_scatter_vacuum_ucoarse_cbc.py @@ -0,0 +1,154 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +transport_2d_cyl_2g_scatter_vacuum_ucoarse.py: 2-group scatter on gmsh coarse +mesh with vacuum boundaries. +Expected: PASS=1 on comparison with orthogonal reference solution computed in +this input. +""" + +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.mesh import FromFileMeshGenerator, OrthogonalMeshGenerator + from pyopensn.logvol import RPPLogicalVolume + from pyopensn.xs import MultiGroupXS + from pyopensn.source import VolumetricSource + from pyopensn.aquad import GLCProductQuadrature2DRZ + from pyopensn.solver import DiscreteOrdinatesCurvilinearProblem, SteadyStateSourceSolver + from pyopensn.fieldfunc import FieldFunctionInterpolationVolume + + +if __name__ == "__main__": + def run_problem( + grid, + xs, + src, + quad, + angle_agg, + inner_method, + l_abs_tol, + l_max_its, + vol, + boundary_conditions, + ): + problem = DiscreteOrdinatesCurvilinearProblem( + mesh=grid, + sweep_type="CBC", + num_groups=2, + groupsets=[ + { + "groups_from_to": (0, 1), + "angular_quadrature": quad, + "angle_aggregation_type": angle_agg, + "inner_linear_method": inner_method, + "l_abs_tol": l_abs_tol, + "l_max_its": l_max_its, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + volumetric_sources=[src], + boundary_conditions=boundary_conditions, + ) + + solver = SteadyStateSourceSolver(problem=problem) + solver.Initialize() + solver.Execute() + + fflist = problem.GetScalarFluxFieldFunction(only_scalar_flux=False) + + ffi_g0 = FieldFunctionInterpolationVolume() + ffi_g0.SetOperationType("max") + ffi_g0.SetLogicalVolume(vol) + ffi_g0.AddFieldFunction(fflist[0][0]) + ffi_g0.Execute() + phi_g0 = ffi_g0.GetValue() + + ffi_g1 = FieldFunctionInterpolationVolume() + ffi_g1.SetOperationType("max") + ffi_g1.SetLogicalVolume(vol) + ffi_g1.AddFieldFunction(fflist[1][0]) + ffi_g1.Execute() + phi_g1 = ffi_g1.GetValue() + + return phi_g0, phi_g1 + + level = "coarse" + mesh_file = "../../../../assets/mesh/rz_tri_coarse.msh" + ref_nr = 200 + ref_nz = 200 + tol = 5.0e-2 + + meshgen = FromFileMeshGenerator(filename=mesh_file, coord_sys="cylindrical") + grid = meshgen.Execute() + vol = RPPLogicalVolume(xmin=0.0, xmax=1.0, ymin=0.0, ymax=2.0, infz=True) + grid.SetUniformBlockID(0) + + xs = MultiGroupXS() + xs.LoadFromOpenSn("../../../../assets/xs/transport_2d_cyl_2g_2.xs") + src = VolumetricSource(block_ids=[0], group_strength=[1.5, 0.8]) + + quad = GLCProductQuadrature2DRZ(n_polar=4, n_azimuthal=8, scattering_order=0) + gmsh_bcs = [ + {"name": "rmin", "type": "reflecting"}, + {"name": "rmax", "type": "vacuum"}, + {"name": "zmin", "type": "vacuum"}, + {"name": "zmax", "type": "vacuum"}, + ] + ref_bcs = [ + {"name": "rmin", "type": "reflecting"}, + {"name": "rmax", "type": "vacuum"}, + {"name": "zmin", "type": "vacuum"}, + {"name": "zmax", "type": "vacuum"}, + ] + phi_g0, phi_g1 = run_problem( + grid, + xs, + src, + quad, + angle_agg="single", + inner_method="petsc_richardson", + l_abs_tol=1.0e-12, + l_max_its=5000, + vol=vol, + boundary_conditions=gmsh_bcs, + ) + + ref_nodes = [ + [i * (1.0 / ref_nr) for i in range(ref_nr + 1)], + [i * (2.0 / ref_nz) for i in range(ref_nz + 1)], + ] + ref_meshgen = OrthogonalMeshGenerator(node_sets=ref_nodes, coord_sys="cylindrical") + ref_grid = ref_meshgen.Execute() + ref_grid.SetBlockIDFromLogicalVolume(vol, 0, True) + phi_ref_g0, phi_ref_g1 = run_problem( + ref_grid, + xs, + src, + quad, + angle_agg="single", + inner_method="petsc_richardson", + l_abs_tol=1.0e-12, + l_max_its=200, + vol=vol, + boundary_conditions=ref_bcs, + ) + + ok_g0 = abs(phi_g0 - phi_ref_g0) / max(abs(phi_ref_g0), 1.0e-12) <= tol + ok_g1 = abs(phi_g1 - phi_ref_g1) / max(abs(phi_ref_g1), 1.0e-12) <= tol + ok = ok_g0 and ok_g1 + + if rank == 0: + print(f"G0_MAX {phi_g0:.12e}") + print(f"G1_MAX {phi_g1:.12e}") + print(f"G0_REF {phi_ref_g0:.12e}") + print(f"G1_REF {phi_ref_g1:.12e}") + print(f"TOL {tol:.6e}") + print(f"PASS {int(ok)}") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_2g_scatter_vacuum_ufine_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_2g_scatter_vacuum_ufine_cbc.py new file mode 100644 index 000000000..2326c2f8e --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_2g_scatter_vacuum_ufine_cbc.py @@ -0,0 +1,154 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +transport_2d_cyl_2g_scatter_vacuum_ufine.py: 2-group scatter on gmsh fine +mesh with vacuum boundaries. +Expected: PASS=1 on comparison with orthogonal reference solution computed in +this input. +""" + +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.mesh import FromFileMeshGenerator, OrthogonalMeshGenerator + from pyopensn.logvol import RPPLogicalVolume + from pyopensn.xs import MultiGroupXS + from pyopensn.source import VolumetricSource + from pyopensn.aquad import GLCProductQuadrature2DRZ + from pyopensn.solver import DiscreteOrdinatesCurvilinearProblem, SteadyStateSourceSolver + from pyopensn.fieldfunc import FieldFunctionInterpolationVolume + + +if __name__ == "__main__": + def run_problem( + grid, + xs, + src, + quad, + angle_agg, + inner_method, + l_abs_tol, + l_max_its, + vol, + boundary_conditions, + ): + problem = DiscreteOrdinatesCurvilinearProblem( + mesh=grid, + sweep_type="CBC", + num_groups=2, + groupsets=[ + { + "groups_from_to": (0, 1), + "angular_quadrature": quad, + "angle_aggregation_type": angle_agg, + "inner_linear_method": inner_method, + "l_abs_tol": l_abs_tol, + "l_max_its": l_max_its, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + volumetric_sources=[src], + boundary_conditions=boundary_conditions, + ) + + solver = SteadyStateSourceSolver(problem=problem) + solver.Initialize() + solver.Execute() + + fflist = problem.GetScalarFluxFieldFunction(only_scalar_flux=False) + + ffi_g0 = FieldFunctionInterpolationVolume() + ffi_g0.SetOperationType("max") + ffi_g0.SetLogicalVolume(vol) + ffi_g0.AddFieldFunction(fflist[0][0]) + ffi_g0.Execute() + phi_g0 = ffi_g0.GetValue() + + ffi_g1 = FieldFunctionInterpolationVolume() + ffi_g1.SetOperationType("max") + ffi_g1.SetLogicalVolume(vol) + ffi_g1.AddFieldFunction(fflist[1][0]) + ffi_g1.Execute() + phi_g1 = ffi_g1.GetValue() + + return phi_g0, phi_g1 + + level = "fine" + mesh_file = "../../../../assets/mesh/rz_tri_fine.msh" + ref_nr = 200 + ref_nz = 200 + tol = 5.0e-2 + + meshgen = FromFileMeshGenerator(filename=mesh_file, coord_sys="cylindrical") + grid = meshgen.Execute() + vol = RPPLogicalVolume(xmin=0.0, xmax=1.0, ymin=0.0, ymax=2.0, infz=True) + grid.SetUniformBlockID(0) + + xs = MultiGroupXS() + xs.LoadFromOpenSn("../../../../assets/xs/transport_2d_cyl_2g_2.xs") + src = VolumetricSource(block_ids=[0], group_strength=[1.5, 0.8]) + + quad = GLCProductQuadrature2DRZ(n_polar=4, n_azimuthal=8, scattering_order=0) + gmsh_bcs = [ + {"name": "rmin", "type": "reflecting"}, + {"name": "rmax", "type": "vacuum"}, + {"name": "zmin", "type": "vacuum"}, + {"name": "zmax", "type": "vacuum"}, + ] + ref_bcs = [ + {"name": "rmin", "type": "reflecting"}, + {"name": "rmax", "type": "vacuum"}, + {"name": "zmin", "type": "vacuum"}, + {"name": "zmax", "type": "vacuum"}, + ] + phi_g0, phi_g1 = run_problem( + grid, + xs, + src, + quad, + angle_agg="single", + inner_method="petsc_richardson", + l_abs_tol=1.0e-12, + l_max_its=5000, + vol=vol, + boundary_conditions=gmsh_bcs, + ) + + ref_nodes = [ + [i * (1.0 / ref_nr) for i in range(ref_nr + 1)], + [i * (2.0 / ref_nz) for i in range(ref_nz + 1)], + ] + ref_meshgen = OrthogonalMeshGenerator(node_sets=ref_nodes, coord_sys="cylindrical") + ref_grid = ref_meshgen.Execute() + ref_grid.SetBlockIDFromLogicalVolume(vol, 0, True) + phi_ref_g0, phi_ref_g1 = run_problem( + ref_grid, + xs, + src, + quad, + angle_agg="single", + inner_method="petsc_richardson", + l_abs_tol=1.0e-12, + l_max_its=200, + vol=vol, + boundary_conditions=ref_bcs, + ) + + ok_g0 = abs(phi_g0 - phi_ref_g0) / max(abs(phi_ref_g0), 1.0e-12) <= tol + ok_g1 = abs(phi_g1 - phi_ref_g1) / max(abs(phi_ref_g1), 1.0e-12) <= tol + ok = ok_g0 and ok_g1 + + if rank == 0: + print(f"G0_MAX {phi_g0:.12e}") + print(f"G1_MAX {phi_g1:.12e}") + print(f"G0_REF {phi_ref_g0:.12e}") + print(f"G1_REF {phi_ref_g1:.12e}") + print(f"TOL {tol:.6e}") + print(f"PASS {int(ok)}") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_2g_scatter_vacuum_umed_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_2g_scatter_vacuum_umed_cbc.py new file mode 100644 index 000000000..062f72062 --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_2g_scatter_vacuum_umed_cbc.py @@ -0,0 +1,154 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +transport_2d_cyl_2g_scatter_vacuum_umed.py: 2-group scatter on gmsh medium +mesh with vacuum boundaries. +Expected: PASS=1 on comparison with orthogonal reference solution computed in +this input. +""" + +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.mesh import FromFileMeshGenerator, OrthogonalMeshGenerator + from pyopensn.logvol import RPPLogicalVolume + from pyopensn.xs import MultiGroupXS + from pyopensn.source import VolumetricSource + from pyopensn.aquad import GLCProductQuadrature2DRZ + from pyopensn.solver import DiscreteOrdinatesCurvilinearProblem, SteadyStateSourceSolver + from pyopensn.fieldfunc import FieldFunctionInterpolationVolume + + +if __name__ == "__main__": + def run_problem( + grid, + xs, + src, + quad, + angle_agg, + inner_method, + l_abs_tol, + l_max_its, + vol, + boundary_conditions, + ): + problem = DiscreteOrdinatesCurvilinearProblem( + mesh=grid, + sweep_type="CBC", + num_groups=2, + groupsets=[ + { + "groups_from_to": (0, 1), + "angular_quadrature": quad, + "angle_aggregation_type": angle_agg, + "inner_linear_method": inner_method, + "l_abs_tol": l_abs_tol, + "l_max_its": l_max_its, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + volumetric_sources=[src], + boundary_conditions=boundary_conditions, + ) + + solver = SteadyStateSourceSolver(problem=problem) + solver.Initialize() + solver.Execute() + + fflist = problem.GetScalarFluxFieldFunction(only_scalar_flux=False) + + ffi_g0 = FieldFunctionInterpolationVolume() + ffi_g0.SetOperationType("max") + ffi_g0.SetLogicalVolume(vol) + ffi_g0.AddFieldFunction(fflist[0][0]) + ffi_g0.Execute() + phi_g0 = ffi_g0.GetValue() + + ffi_g1 = FieldFunctionInterpolationVolume() + ffi_g1.SetOperationType("max") + ffi_g1.SetLogicalVolume(vol) + ffi_g1.AddFieldFunction(fflist[1][0]) + ffi_g1.Execute() + phi_g1 = ffi_g1.GetValue() + + return phi_g0, phi_g1 + + level = "medium" + mesh_file = "../../../../assets/mesh/rz_tri_medium.msh" + ref_nr = 200 + ref_nz = 200 + tol = 5.0e-2 + + meshgen = FromFileMeshGenerator(filename=mesh_file, coord_sys="cylindrical") + grid = meshgen.Execute() + vol = RPPLogicalVolume(xmin=0.0, xmax=1.0, ymin=0.0, ymax=2.0, infz=True) + grid.SetUniformBlockID(0) + + xs = MultiGroupXS() + xs.LoadFromOpenSn("../../../../assets/xs/transport_2d_cyl_2g_2.xs") + src = VolumetricSource(block_ids=[0], group_strength=[1.5, 0.8]) + + quad = GLCProductQuadrature2DRZ(n_polar=4, n_azimuthal=8, scattering_order=0) + gmsh_bcs = [ + {"name": "rmin", "type": "reflecting"}, + {"name": "rmax", "type": "vacuum"}, + {"name": "zmin", "type": "vacuum"}, + {"name": "zmax", "type": "vacuum"}, + ] + ref_bcs = [ + {"name": "rmin", "type": "reflecting"}, + {"name": "rmax", "type": "vacuum"}, + {"name": "zmin", "type": "vacuum"}, + {"name": "zmax", "type": "vacuum"}, + ] + phi_g0, phi_g1 = run_problem( + grid, + xs, + src, + quad, + angle_agg="single", + inner_method="petsc_richardson", + l_abs_tol=1.0e-12, + l_max_its=5000, + vol=vol, + boundary_conditions=gmsh_bcs, + ) + + ref_nodes = [ + [i * (1.0 / ref_nr) for i in range(ref_nr + 1)], + [i * (2.0 / ref_nz) for i in range(ref_nz + 1)], + ] + ref_meshgen = OrthogonalMeshGenerator(node_sets=ref_nodes, coord_sys="cylindrical") + ref_grid = ref_meshgen.Execute() + ref_grid.SetBlockIDFromLogicalVolume(vol, 0, True) + phi_ref_g0, phi_ref_g1 = run_problem( + ref_grid, + xs, + src, + quad, + angle_agg="single", + inner_method="petsc_richardson", + l_abs_tol=1.0e-12, + l_max_its=200, + vol=vol, + boundary_conditions=ref_bcs, + ) + + ok_g0 = abs(phi_g0 - phi_ref_g0) / max(abs(phi_ref_g0), 1.0e-12) <= tol + ok_g1 = abs(phi_g1 - phi_ref_g1) / max(abs(phi_ref_g1), 1.0e-12) <= tol + ok = ok_g0 and ok_g1 + + if rank == 0: + print(f"G0_MAX {phi_g0:.12e}") + print(f"G1_MAX {phi_g1:.12e}") + print(f"G0_REF {phi_ref_g0:.12e}") + print(f"G1_REF {phi_ref_g1:.12e}") + print(f"TOL {tol:.6e}") + print(f"PASS {int(ok)}") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_annulus_zero_solution_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_annulus_zero_solution_cbc.py new file mode 100644 index 000000000..a39136e4d --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_annulus_zero_solution_cbc.py @@ -0,0 +1,72 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +transport_2d_cyl_annulus_zero_solution.py: Annulus with vacuum and zero source +Expected: MAX=0.0 from analytic zero-source solution +""" + +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.mesh import FromFileMeshGenerator + from pyopensn.logvol import RPPLogicalVolume + from pyopensn.xs import MultiGroupXS + from pyopensn.aquad import GLCProductQuadrature2DRZ + from pyopensn.solver import DiscreteOrdinatesCurvilinearProblem, SteadyStateSourceSolver + from pyopensn.fieldfunc import FieldFunctionInterpolationVolume + + +if __name__ == "__main__": + if size != 4: + sys.exit(f"Incorrect number of processors. Expected 4 but got {size}.") + + meshgen = FromFileMeshGenerator( + filename="../../../../assets/mesh/rz_annulus_single.msh", coord_sys="cylindrical" + ) + grid = meshgen.Execute() + + vol = RPPLogicalVolume(xmin=0.25, xmax=1.0, ymin=0.0, ymax=2.0, infz=True) + + xs = MultiGroupXS() + xs.CreateSimpleOneGroup(1.0, 0.0) + + quad = GLCProductQuadrature2DRZ(n_polar=4, n_azimuthal=8, scattering_order=0) + problem = DiscreteOrdinatesCurvilinearProblem( + mesh=grid, + sweep_type="CBC", + num_groups=1, + groupsets=[ + { + "groups_from_to": (0, 0), + "angular_quadrature": quad, + "angle_aggregation_type": "azimuthal", + "inner_linear_method": "petsc_gmres", + "l_abs_tol": 1.0e-12, + "l_max_its": 100, + } + ], + xs_map=[{"block_ids": [1], "xs": xs}], + ) + + solver = SteadyStateSourceSolver(problem=problem) + solver.Initialize() + solver.Execute() + + ff = problem.GetScalarFluxFieldFunction(only_scalar_flux=False)[0][0] + + ffi_max = FieldFunctionInterpolationVolume() + ffi_max.SetOperationType("max") + ffi_max.SetLogicalVolume(vol) + ffi_max.AddFieldFunction(ff) + ffi_max.Execute() + phi_max = ffi_max.GetValue() + + if rank == 0: + print(f"MAX {phi_max:.12e}") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_disk_leakage_reflecting_bcs_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_disk_leakage_reflecting_bcs_cbc.py new file mode 100644 index 000000000..c7c3c50f2 --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_disk_leakage_reflecting_bcs_cbc.py @@ -0,0 +1,103 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +transport_2d_cyl_disk_leakage_reflecting_bcs.py: Disk with radial vacuum +boundary and reflecting z boundaries. +Expected: From 3D Cartesian disk problem +(transport_3d_disk_leakage_reflecting_bcs.py): +Top=-6.190966583965322e-15 +Bottom=-5.231802747012948e-16 +Side=97.24938933424164 +Total=97.24938933424164 +""" + +import os +import sys +import math + +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.mesh import OrthogonalMeshGenerator + from pyopensn.xs import MultiGroupXS + from pyopensn.source import VolumetricSource + from pyopensn.aquad import GLCProductQuadrature2DRZ + from pyopensn.solver import DiscreteOrdinatesCurvilinearProblem, SteadyStateSourceSolver + + +if __name__ == "__main__": + rmin = 0.0 + rmax = 0.2 + zmin = 0.0 + zmax = 0.1 + nr = 100 + nz = 50 + r_nodes = [rmin + i * ((rmax - rmin) / nr) for i in range(nr + 1)] + z_nodes = [zmin + i * ((zmax - zmin) / nz) for i in range(nz + 1)] + meshgen = OrthogonalMeshGenerator(node_sets=[r_nodes, z_nodes], coord_sys="cylindrical") + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + + xs = MultiGroupXS() + xs.CreateSimpleOneGroup(sigma_t=1.0, c=0.0) + + vol_src = VolumetricSource(block_ids=[0], group_strength=[9754.5]) + + quad = GLCProductQuadrature2DRZ(n_polar=16, n_azimuthal=32, scattering_order=0) + bcs = [ + {"name": "rmin", "type": "reflecting"}, # axis symmetry + {"name": "rmax", "type": "vacuum"}, + {"name": "zmin", "type": "reflecting"}, + {"name": "zmax", "type": "reflecting"}, + ] + + phys = DiscreteOrdinatesCurvilinearProblem( + mesh=grid, + sweep_type="CBC", + num_groups=1, + groupsets=[ + { + "groups_from_to": (0, 0), + "angular_quadrature": quad, + "angle_aggregation_type": "single", + "inner_linear_method": "petsc_gmres", + "l_abs_tol": 1.0e-6, + "l_max_its": 50, + "gmres_restart_interval": 10, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + boundary_conditions=bcs, + volumetric_sources=[vol_src], + options={"save_angular_flux": True}, + ) + + solver = SteadyStateSourceSolver(problem=phys) + solver.Initialize() + solver.Execute() + + leakage = phys.ComputeLeakage(["rmax", "zmin", "zmax"]) + lkg_side = leakage["rmax"] + lkg_bottom = leakage["zmin"] + lkg_top = leakage["zmax"] + lkg_total = lkg_side + lkg_top + lkg_bottom + scale = 2.0 * math.pi + + top_val = float(lkg_top.item()) * scale + bottom_val = float(lkg_bottom.item()) * scale + side_val = float(lkg_side.item()) * scale + total_val = float(lkg_total.item()) * scale + + if rank == 0: + volume = math.pi * rmax * rmax * (zmax - zmin) + total_source = 9754.5 * volume + print(f"Top leakage={top_val}") + print(f"Bottom leakage={bottom_val}") + print(f"Side leakage={side_val}") + print(f"Total leakage={total_val}") + print(f"Total source={total_source}") + if total_source > 0.0: + print(f"Leakage/Source={total_val / total_source}") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_disk_leakage_vacuum_bcs_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_disk_leakage_vacuum_bcs_cbc.py new file mode 100644 index 000000000..3aaa8bdb5 --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_disk_leakage_vacuum_bcs_cbc.py @@ -0,0 +1,103 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +transport_2d_cyl_disk_leakage_vacuum_bcs.py: Disk with all vacuum boundaries +Expected: From 3D Cartesian disk problem (transport_3d_disk_leakage_vacuum_bcs.py): +Top=36.173790857624596 +Bottom=36.17379086382096 +Side=39.74878117223803 +Total=112.0963628936836 +""" + +import os +import sys +import math + +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.mesh import OrthogonalMeshGenerator + from pyopensn.xs import MultiGroupXS + from pyopensn.source import VolumetricSource + from pyopensn.aquad import GLCProductQuadrature2DRZ + from pyopensn.solver import DiscreteOrdinatesCurvilinearProblem, SteadyStateSourceSolver + + +if __name__ == "__main__": + rmin = 0.0 + rmax = 0.2 + zmin = 0.0 + zmax = 0.1 + nr = 100 + nz = 50 + r_nodes = [rmin + i * ((rmax - rmin) / nr) for i in range(nr + 1)] + z_nodes = [zmin + i * ((zmax - zmin) / nz) for i in range(nz + 1)] + meshgen = OrthogonalMeshGenerator(node_sets=[r_nodes, z_nodes], coord_sys="cylindrical") + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + + xs = MultiGroupXS() + xs.CreateSimpleOneGroup(sigma_t=1.0, c=0.0) + + # Match 3D test source: Q=122.58 particles/s over V=pi*r^2*h + # Strength per unit volume: 122.58 / (pi*0.2^2*0.1) = 9754.5 + vol_src = VolumetricSource(block_ids=[0], group_strength=[9754.5]) + + quad = GLCProductQuadrature2DRZ(n_polar=16, n_azimuthal=32, scattering_order=0) + bcs = [ + {"name": "rmin", "type": "reflecting"}, # axis symmetry + {"name": "rmax", "type": "vacuum"}, + {"name": "zmin", "type": "vacuum"}, + {"name": "zmax", "type": "vacuum"}, + ] + + phys = DiscreteOrdinatesCurvilinearProblem( + mesh=grid, + sweep_type="CBC", + num_groups=1, + groupsets=[ + { + "groups_from_to": (0, 0), + "angular_quadrature": quad, + "angle_aggregation_type": "single", + "inner_linear_method": "petsc_gmres", + "l_abs_tol": 1.0e-6, + "l_max_its": 30, + "gmres_restart_interval": 10, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + boundary_conditions=bcs, + volumetric_sources=[vol_src], + options={"save_angular_flux": True}, + ) + + solver = SteadyStateSourceSolver(problem=phys) + solver.Initialize() + solver.Execute() + + leakage = phys.ComputeLeakage(["rmax", "zmin", "zmax"]) + lkg_side = leakage["rmax"] + lkg_bottom = leakage["zmin"] + lkg_top = leakage["zmax"] + lkg_total = lkg_side + lkg_top + lkg_bottom + scale = 2.0 * math.pi + + top_val = float(lkg_top.item()) * scale + bottom_val = float(lkg_bottom.item()) * scale + side_val = float(lkg_side.item()) * scale + total_val = float(lkg_total.item()) * scale + + if rank == 0: + volume = math.pi * rmax * rmax * (zmax - zmin) + total_source = 9754.5 * volume + print(f"Top leakage={top_val}") + print(f"Bottom leakage={bottom_val}") + print(f"Side leakage={side_val}") + print(f"Total leakage={total_val}") + print(f"Total source={total_source}") + if total_source > 0.0: + print(f"Leakage/Source={total_val / total_source}") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_disk_source_vacuum_hdpe_struct_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_disk_source_vacuum_hdpe_struct_cbc.py new file mode 100644 index 000000000..f5adba9b5 --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_disk_source_vacuum_hdpe_struct_cbc.py @@ -0,0 +1,86 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +transport_2d_cyl_disk_source_vacuum_hdpe_struct.py: Disk with HDPE OpenMC cross sections +Expected: PHI_MAX=1.224546944901e-01 from 3D Cartesian reference problem +""" + +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.mesh import OrthogonalMeshGenerator + from pyopensn.logvol import RPPLogicalVolume + from pyopensn.xs import MultiGroupXS + from pyopensn.source import VolumetricSource + from pyopensn.aquad import GLCProductQuadrature2DRZ + from pyopensn.solver import DiscreteOrdinatesCurvilinearProblem, SteadyStateSourceSolver + from pyopensn.fieldfunc import FieldFunctionInterpolationVolume + + +if __name__ == "__main__": + # Structured cylindrical mesh on 0<=r<=0.2, 0<=z<=0.1. + nr = 100 + nz = 50 + r_nodes = [0.2 * i / nr for i in range(nr + 1)] + z_nodes = [0.1 * i / nz for i in range(nz + 1)] + meshgen = OrthogonalMeshGenerator(node_sets=[r_nodes, z_nodes], coord_sys="cylindrical") + grid = meshgen.Execute() + + vol = RPPLogicalVolume(xmin=0.0, xmax=0.2, ymin=0.0, ymax=0.1, infz=True) + grid.SetBlockIDFromLogicalVolume(vol, 0, True) + + xs = MultiGroupXS() + xs.LoadFromOpenMC("../../../../assets/xs/HDPE.h5", "set1", 294.0) + + num_groups = 172 + strength = [0.0 for _ in range(num_groups)] + strength[0] = 1.0 + src = VolumetricSource(block_ids=[0], group_strength=strength) + + quad = GLCProductQuadrature2DRZ(n_polar=8, n_azimuthal=16, scattering_order=0) + problem = DiscreteOrdinatesCurvilinearProblem( + mesh=grid, + sweep_type="CBC", + num_groups=num_groups, + groupsets=[ + { + "groups_from_to": (0, num_groups - 1), + "angular_quadrature": quad, + "angle_aggregation_type": "azimuthal", + "inner_linear_method": "petsc_gmres", + "l_abs_tol": 1.0e-8, + "l_max_its": 300, + "gmres_restart_interval": 100, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + volumetric_sources=[src], + boundary_conditions=[ + {"name": "rmin", "type": "reflecting"}, + {"name": "rmax", "type": "vacuum"}, + {"name": "zmin", "type": "vacuum"}, + {"name": "zmax", "type": "vacuum"}, + ], + ) + + solver = SteadyStateSourceSolver(problem=problem) + solver.Initialize() + solver.Execute() + + fflist = problem.GetScalarFluxFieldFunction(only_scalar_flux=False) + ffi = FieldFunctionInterpolationVolume() + ffi.SetOperationType("max") + ffi.SetLogicalVolume(vol) + ffi.AddFieldFunction(fflist[0][0]) + ffi.Execute() + phi_max = ffi.GetValue() + + if rank == 0: + print(f"PHI_MAX {phi_max:.12e}") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_reflective_uniform_1g_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_reflective_uniform_1g_cbc.py new file mode 100644 index 000000000..7dddac425 --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_reflective_uniform_1g_cbc.py @@ -0,0 +1,86 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +transport_2d_cyl_reflective_uniform_1g.py: Uniform 1-group RZ with rmin reflecting +Expected: MAX=0.8437473489223 from 3D Cartesian reference problem +""" + +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.mesh import FromFileMeshGenerator + from pyopensn.logvol import RPPLogicalVolume + from pyopensn.xs import MultiGroupXS + from pyopensn.source import VolumetricSource + from pyopensn.aquad import GLCProductQuadrature2DRZ + from pyopensn.solver import DiscreteOrdinatesCurvilinearProblem, SteadyStateSourceSolver + from pyopensn.fieldfunc import FieldFunctionInterpolationVolume + + +if __name__ == "__main__": + if size != 4: + sys.exit(f"Incorrect number of processors. Expected 4 but got {size}.") + + meshgen = FromFileMeshGenerator( + filename="../../../../assets/mesh/rz_rect_single.msh", coord_sys="cylindrical" + ) + grid = meshgen.Execute() + + vol = RPPLogicalVolume(xmin=0.0, xmax=1.0, ymin=0.0, ymax=2.0, infz=True) + grid.SetBlockIDFromLogicalVolume(vol, 0, True) + + xs = MultiGroupXS() + sigma_t = 2.0 + c = 0.25 + xs.CreateSimpleOneGroup(sigma_t, c) + + sigma_a = sigma_t * (1.0 - c) + src = VolumetricSource(block_ids=[0], group_strength=[sigma_a]) + + quad = GLCProductQuadrature2DRZ(n_polar=4, n_azimuthal=8, scattering_order=0) + problem = DiscreteOrdinatesCurvilinearProblem( + mesh=grid, + sweep_type="CBC", + num_groups=1, + groupsets=[ + { + "groups_from_to": (0, 0), + "angular_quadrature": quad, + "angle_aggregation_type": "single", + "inner_linear_method": "petsc_richardson", + "l_abs_tol": 1.0e-12, + "l_max_its": 100, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + volumetric_sources=[src], + boundary_conditions=[ + {"name": "rmin", "type": "reflecting"}, + {"name": "rmax", "type": "vacuum"}, + {"name": "zmin", "type": "vacuum"}, + {"name": "zmax", "type": "vacuum"}, + ], + ) + + solver = SteadyStateSourceSolver(problem=problem) + solver.Initialize() + solver.Execute() + + ff = problem.GetScalarFluxFieldFunction(only_scalar_flux=False)[0][0] + + ffi_max = FieldFunctionInterpolationVolume() + ffi_max.SetOperationType("max") + ffi_max.SetLogicalVolume(vol) + ffi_max.AddFieldFunction(ff) + ffi_max.Execute() + phi_max = ffi_max.GetValue() + + if rank == 0: + print(f"MAX {phi_max:.12e}") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_reflective_uniform_2g_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_reflective_uniform_2g_cbc.py new file mode 100644 index 000000000..d395e9cb2 --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_reflective_uniform_2g_cbc.py @@ -0,0 +1,93 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +transport_2d_cyl_reflective_uniform_2g.py: Uniform 2-group with rmin reflecting. +Expected: G0_MAX=0.8437473489224, G1_MAX=0.3658623519312 from 3D Cartesian reference problem +""" + +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.mesh import FromFileMeshGenerator + from pyopensn.logvol import RPPLogicalVolume + from pyopensn.xs import MultiGroupXS + from pyopensn.source import VolumetricSource + from pyopensn.aquad import GLCProductQuadrature2DRZ + from pyopensn.solver import DiscreteOrdinatesCurvilinearProblem, SteadyStateSourceSolver + from pyopensn.fieldfunc import FieldFunctionInterpolationVolume + + +if __name__ == "__main__": + if size != 4: + sys.exit(f"Incorrect number of processors. Expected 4 but got {size}.") + + meshgen = FromFileMeshGenerator( + filename="../../../../assets/mesh/rz_rect_single.msh", coord_sys="cylindrical" + ) + grid = meshgen.Execute() + + vol = RPPLogicalVolume(xmin=0.0, xmax=1.0, ymin=0.0, ymax=2.0, infz=True) + grid.SetBlockIDFromLogicalVolume(vol, 0, True) + + xs = MultiGroupXS() + xs.LoadFromOpenSn("../../../../assets/xs/transport_2d_cyl_2g_2.xs") + + # q0 = sigma_a0 * phi0 = 1.5 * 1.0 + # q1 = sigma_a1 * phi1 - sigma_s(0->1) * phi0 = 2.5 * 0.4 - 0.2 * 1.0 + src = VolumetricSource(block_ids=[0], group_strength=[1.5, 0.8]) + + quad = GLCProductQuadrature2DRZ(n_polar=4, n_azimuthal=8, scattering_order=0) + problem = DiscreteOrdinatesCurvilinearProblem( + mesh=grid, + sweep_type="CBC", + num_groups=2, + groupsets=[ + { + "groups_from_to": (0, 1), + "angular_quadrature": quad, + "angle_aggregation_type": "single", + "inner_linear_method": "petsc_gmres", + "l_abs_tol": 1.0e-12, + "l_max_its": 400, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + volumetric_sources=[src], + boundary_conditions=[ + {"name": "rmin", "type": "reflecting"}, + {"name": "rmax", "type": "vacuum"}, + {"name": "zmin", "type": "vacuum"}, + {"name": "zmax", "type": "vacuum"}, + ], + ) + + solver = SteadyStateSourceSolver(problem=problem) + solver.Initialize() + solver.Execute() + + fflist = problem.GetScalarFluxFieldFunction(only_scalar_flux=False) + + ffi_g0 = FieldFunctionInterpolationVolume() + ffi_g0.SetOperationType("max") + ffi_g0.SetLogicalVolume(vol) + ffi_g0.AddFieldFunction(fflist[0][0]) + ffi_g0.Execute() + phi_g0 = ffi_g0.GetValue() + + ffi_g1 = FieldFunctionInterpolationVolume() + ffi_g1.SetOperationType("max") + ffi_g1.SetLogicalVolume(vol) + ffi_g1.AddFieldFunction(fflist[1][0]) + ffi_g1.Execute() + phi_g1 = ffi_g1.GetValue() + + if rank == 0: + print(f"G0_MAX {phi_g0:.12e}") + print(f"G1_MAX {phi_g1:.12e}") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_unstructured_reflective_uniform_1g_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_unstructured_reflective_uniform_1g_cbc.py new file mode 100644 index 000000000..1d2cf951d --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_unstructured_reflective_uniform_1g_cbc.py @@ -0,0 +1,87 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +transport_2d_cyl_unstructured_reflective_uniform_1g.py: Unstructured 1-group with vacuum boundaries +Expected: MAX=0.9065113759285 from 3D Cartesian reference problem +""" + +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.mesh import FromFileMeshGenerator + from pyopensn.logvol import RPPLogicalVolume + from pyopensn.xs import MultiGroupXS + from pyopensn.source import VolumetricSource + from pyopensn.aquad import GLCProductQuadrature2DRZ + from pyopensn.solver import DiscreteOrdinatesCurvilinearProblem, SteadyStateSourceSolver + from pyopensn.fieldfunc import FieldFunctionInterpolationVolume + + +if __name__ == "__main__": + if size != 4: + sys.exit(f"Incorrect number of processors. Expected 4 but got {size}.") + + meshgen = FromFileMeshGenerator( + filename="../../../../assets/mesh/rz_rect_unstructured.msh", + coord_sys="cylindrical", + ) + grid = meshgen.Execute() + + vol = RPPLogicalVolume(xmin=0.0, xmax=1.0, ymin=0.0, ymax=2.0, infz=True) + grid.SetBlockIDFromLogicalVolume(vol, 0, True) + + xs = MultiGroupXS() + sigma_t = 5.0 + c = 0.8 + xs.CreateSimpleOneGroup(sigma_t, c) + + sigma_a = sigma_t * (1.0 - c) + src = VolumetricSource(block_ids=[0], group_strength=[sigma_a]) + + quad = GLCProductQuadrature2DRZ(n_polar=6, n_azimuthal=12, scattering_order=0) + problem = DiscreteOrdinatesCurvilinearProblem( + mesh=grid, + sweep_type="CBC", + num_groups=1, + groupsets=[ + { + "groups_from_to": (0, 0), + "angular_quadrature": quad, + "angle_aggregation_type": "single", + "inner_linear_method": "petsc_richardson", + "l_abs_tol": 1.0e-12, + "l_max_its": 100, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + volumetric_sources=[src], + boundary_conditions=[ + {"name": "rmin", "type": "reflecting"}, + {"name": "rmax", "type": "vacuum"}, + {"name": "zmin", "type": "vacuum"}, + {"name": "zmax", "type": "vacuum"}, + ], + ) + + solver = SteadyStateSourceSolver(problem=problem) + solver.Initialize() + solver.Execute() + + ff = problem.GetScalarFluxFieldFunction(only_scalar_flux=False)[0][0] + + ffi_max = FieldFunctionInterpolationVolume() + ffi_max.SetOperationType("max") + ffi_max.SetLogicalVolume(vol) + ffi_max.AddFieldFunction(ff) + ffi_max.Execute() + phi_max = ffi_max.GetValue() + + if rank == 0: + print(f"MAX {phi_max:.12e}") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_zero_source_vacuum_1g_struct_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_zero_source_vacuum_1g_struct_cbc.py new file mode 100644 index 000000000..0dc033e0f --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_steady_cyl/transport_2d_cyl_zero_source_vacuum_1g_struct_cbc.py @@ -0,0 +1,81 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +transport_2d_cyl_zero_source_vacuum_1g_struct.py: Zero source, vacuum boundaries +Expected: PHI_MAX=0.0 from analytic zero-source solution +""" + +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.mesh import OrthogonalMeshGenerator + from pyopensn.logvol import RPPLogicalVolume + from pyopensn.xs import MultiGroupXS + from pyopensn.source import VolumetricSource + from pyopensn.aquad import GLCProductQuadrature2DRZ + from pyopensn.solver import DiscreteOrdinatesCurvilinearProblem, SteadyStateSourceSolver + from pyopensn.fieldfunc import FieldFunctionInterpolationVolume + + +if __name__ == "__main__": + nr = 40 + nz = 80 + r_nodes = [i * (1.0 / nr) for i in range(nr + 1)] + z_nodes = [i * (2.0 / nz) for i in range(nz + 1)] + meshgen = OrthogonalMeshGenerator(node_sets=[r_nodes, z_nodes], coord_sys="cylindrical") + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + + xs = MultiGroupXS() + xs.CreateSimpleOneGroup(sigma_t=1.0, c=0.0) + src = VolumetricSource(block_ids=[0], group_strength=[0.0]) + + quad = GLCProductQuadrature2DRZ(n_polar=4, n_azimuthal=8, scattering_order=0) + bcs = [ + {"name": "rmin", "type": "vacuum"}, + {"name": "rmax", "type": "vacuum"}, + {"name": "zmin", "type": "vacuum"}, + {"name": "zmax", "type": "vacuum"}, + ] + + problem = DiscreteOrdinatesCurvilinearProblem( + mesh=grid, + sweep_type="CBC", + num_groups=1, + groupsets=[ + { + "groups_from_to": (0, 0), + "angular_quadrature": quad, + "angle_aggregation_type": "azimuthal", + "inner_linear_method": "petsc_richardson", + "l_abs_tol": 1.0e-12, + "l_max_its": 50, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + volumetric_sources=[src], + boundary_conditions=bcs, + ) + + solver = SteadyStateSourceSolver(problem=problem) + solver.Initialize() + solver.Execute() + + vol = RPPLogicalVolume(xmin=0.0, xmax=1.0, ymin=0.0, ymax=2.0, infz=True) + fflist = problem.GetScalarFluxFieldFunction(only_scalar_flux=False) + + ffi_max = FieldFunctionInterpolationVolume() + ffi_max.SetOperationType("max") + ffi_max.SetLogicalVolume(vol) + ffi_max.AddFieldFunction(fflist[0][0]) + ffi_max.Execute() + phi_max = ffi_max.GetValue() + + if rank == 0: + print(f"PHI_MAX {phi_max:.12e}")