From 9297acf0a61167a04d5bce145809bd021cdb22ac Mon Sep 17 00:00:00 2001 From: "W. Daryl Hawkins" Date: Sun, 18 May 2025 16:49:06 -0500 Subject: [PATCH 1/2] Adding SLEPc solver for k-eigenvalue problems. --- CMakeLists.txt | 12 +- cmake/FindPETSc.cmake | 317 +++++++++++++-- .../math/linear_solver/linear_eigen_solver.h | 11 +- .../slepc_linear_eigen_solver.cc | 103 +++++ .../linear_solver/slepc_linear_eigen_solver.h | 82 ++++ .../slepc_linear_keigen_solver.cc | 373 ++++++++++++++++++ .../slepc_linear_keigen_solver.h | 52 +++ .../solvers/slepc_keigen_solver.cc | 89 +++++ .../solvers/slepc_keigen_solver.h | 36 ++ python/lib/py_app.cc | 1 + python/lib/py_wrappers.h | 1 + python/lib/solver.cc | 44 +++ python/main.cc | 56 ++- ...envalue_transport_3d_5g_upscatter_slepc.py | 236 +++++++++++ .../keigenvalue_transport_3d_slepc_U235.py | 98 +++++ ...igenvalue_transport_3d_slepc_power_U235.py | 98 +++++ .../simple_5g_upscatter_fissile.xs | 45 ++- .../transport_keigen/tests.json | 66 +++- 18 files changed, 1649 insertions(+), 71 deletions(-) create mode 100644 framework/math/linear_solver/slepc_linear_eigen_solver.cc create mode 100644 framework/math/linear_solver/slepc_linear_eigen_solver.h create mode 100644 modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/slepc_linear_keigen_solver.cc create mode 100644 modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/slepc_linear_keigen_solver.h create mode 100644 modules/linear_boltzmann_solvers/discrete_ordinates_problem/solvers/slepc_keigen_solver.cc create mode 100644 modules/linear_boltzmann_solvers/discrete_ordinates_problem/solvers/slepc_keigen_solver.h create mode 100644 test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_5g_upscatter_slepc.py create mode 100644 test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_slepc_U235.py create mode 100644 test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_slepc_power_U235.py diff --git a/CMakeLists.txt b/CMakeLists.txt index efd3b2037a..3c74d61fb6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -191,14 +191,14 @@ endif() find_package(PETSc 3.17 REQUIRED) # Check PETSc uses 64bit indices -check_symbol_exists(PETSC_USE_64BIT_INDICES "${PETSC_INCLUDE_DIR}/petscconf.h" PETSC_USE_64BIT_INDICES) -if(NOT ${PETSC_USE_64BIT_INDICES} MATCHES 1) +check_symbol_exists(PETSC_USE_64BIT_INDICES "${PETSC_CONF_HEADER}" PETSC_USE_64BIT_INDICES) +if(NOT PETSC_USE_64BIT_INDICES) message(FATAL_ERROR "PETSc has not been configured with the flag --with-64-bit-indices\n") endif() # Check PETSc includes support for PTSCOTCH -check_symbol_exists(PETSC_HAVE_PTSCOTCH "${PETSC_INCLUDE_DIR}/petscconf.h" PETSC_HAVE_PTSCOTCH) -if(NOT ${PETSC_HAVE_PTSCOTCH} MATCHES 1) +check_symbol_exists(PETSC_HAVE_PTSCOTCH "${PETSC_CONF_HEADER}" PETSC_HAVE_PTSCOTCH) +if(NOT PETSC_HAVE_PTSCOTCH) message(FATAL_ERROR "PETSc has not been configured with PTSCOTCH\n") endif() @@ -209,13 +209,13 @@ if(NOT ${PETSC_HAVE_SUPERLU_DIST} MATCHES 1 OR NOT PETSC_HAVE_SUPERLU_DIST) endif() # Check PetscScalar is double -check_symbol_exists(PETSC_USE_REAL_DOUBLE "${PETSC_INCLUDE_DIR}/petscconf.h" PETSC_USE_REAL_DOUBLE) +check_symbol_exists(PETSC_USE_REAL_DOUBLE "${PETSC_CONF_HEADER}" PETSC_USE_REAL_DOUBLE) if(NOT PETSC_USE_REAL_DOUBLE) message(FATAL_ERROR "PETSc has not been configured with double precision reals.\n" "Reconfigure PETSc with --with-precision=double") endif() -check_symbol_exists(PETSC_USE_COMPLEX "${PETSC_INCLUDE_DIR}/petscconf.h" PETSC_USE_COMPLEX) +check_symbol_exists(PETSC_USE_COMPLEX "${PETSC_CONF_HEADER}" PETSC_USE_COMPLEX) if(PETSC_USE_COMPLEX) message(FATAL_ERROR "PETSc has been configured with complex scalars.\n" diff --git a/cmake/FindPETSc.cmake b/cmake/FindPETSc.cmake index e3e1668787..13b824e261 100644 --- a/cmake/FindPETSc.cmake +++ b/cmake/FindPETSc.cmake @@ -1,50 +1,293 @@ -# Find PETSc +# FindPETSc.cmake # # Once done this will define -# PETSC_FOUND - System has PETSc -# PETSC_INCLUDE_DIR - The PETSc include directory -# PETSC_LIBRARY - The PETSc library +# PETSc_FOUND - System has PETSc and SLEPc +# PETSC_FOUND - Compatibility alias for PETSc_FOUND +# PETSC_INCLUDE_DIR - PETSc/SLEPc include directories +# PETSC_LIBRARY - PETSc/SLEPc libraries or imported targets # PETSC_VERSION - The PETSc version +# PETSC_CONF_HEADER - Full path to petscconf.h +# SLEPC_FOUND - System has SLEPc +# SLEPc_FOUND - Compatibility alias for SLEPC_FOUND +# SLEPC_INCLUDE_DIR - The SLEPc include directory +# SLEPC_LIBRARY - The SLEPc library or imported target +# SLEPC_VERSION - The SLEPc version -find_path( - PETSC_INCLUDE_DIR - petsc.h - PATHS - $ENV{PETSC_DIR}/include - $ENV{PETSC_ROOT}/include +include(FindPackageHandleStandardArgs) + +set(_petsc_roots) +foreach(_var IN ITEMS PETSC_DIR PETSC_ROOT) + if(DEFINED ENV{${_var}} AND NOT "$ENV{${_var}}" STREQUAL "") + list(APPEND _petsc_roots "$ENV{${_var}}") + endif() +endforeach() +list(REMOVE_DUPLICATES _petsc_roots) + +set(_slepc_roots) +foreach(_var IN ITEMS SLEPC_DIR SLEPC_ROOT) + if(DEFINED ENV{${_var}} AND NOT "$ENV{${_var}}" STREQUAL "") + list(APPEND _slepc_roots "$ENV{${_var}}") + endif() +endforeach() +list(REMOVE_DUPLICATES _slepc_roots) + +set(_petsc_config_paths) +set(_petsc_include_hints) +set(_petsc_library_hints) +foreach(_root IN LISTS _petsc_roots) + list(APPEND _petsc_config_paths "${_root}/lib/cmake/PETSc") + list(APPEND _petsc_include_hints "${_root}/include") + list(APPEND _petsc_library_hints "${_root}/lib") + if(DEFINED ENV{PETSC_ARCH} AND NOT "$ENV{PETSC_ARCH}" STREQUAL "") + list(APPEND _petsc_config_paths "${_root}/$ENV{PETSC_ARCH}/lib/cmake/PETSc") + list(APPEND _petsc_include_hints "${_root}/$ENV{PETSC_ARCH}/include") + list(APPEND _petsc_library_hints "${_root}/$ENV{PETSC_ARCH}/lib") + endif() +endforeach() + +if(DEFINED PETSC_INCLUDE_DIR AND NOT "${PETSC_INCLUDE_DIR}" MATCHES "-NOTFOUND$") + list(APPEND _petsc_include_hints ${PETSC_INCLUDE_DIR}) +endif() +if(DEFINED PETSC_LIBRARY AND NOT "${PETSC_LIBRARY}" MATCHES "-NOTFOUND$") + foreach(_library IN LISTS PETSC_LIBRARY) + if(EXISTS "${_library}") + get_filename_component(_library_dir "${_library}" DIRECTORY) + get_filename_component(_prefix_dir "${_library_dir}" DIRECTORY) + list(APPEND _petsc_library_hints "${_library_dir}") + list(APPEND _petsc_include_hints "${_prefix_dir}/include") + endif() + endforeach() +endif() +list(REMOVE_DUPLICATES _petsc_include_hints) +list(REMOVE_DUPLICATES _petsc_library_hints) + +set(_slepc_config_paths) +set(_slepc_include_hints) +set(_slepc_library_hints) +foreach(_root IN LISTS _slepc_roots) + list(APPEND _slepc_config_paths "${_root}/lib/cmake/SLEPc") + list(APPEND _slepc_include_hints "${_root}/include") + list(APPEND _slepc_library_hints "${_root}/lib") + if(DEFINED ENV{SLEPC_ARCH} AND NOT "$ENV{SLEPC_ARCH}" STREQUAL "") + list(APPEND _slepc_config_paths "${_root}/$ENV{SLEPC_ARCH}/lib/cmake/SLEPc") + list(APPEND _slepc_include_hints "${_root}/$ENV{SLEPC_ARCH}/include") + list(APPEND _slepc_library_hints "${_root}/$ENV{SLEPC_ARCH}/lib") + endif() +endforeach() +foreach(_root IN LISTS _petsc_roots) + list(APPEND _slepc_include_hints "${_root}/include/slepc") + list(APPEND _slepc_library_hints "${_root}/lib") + if(DEFINED ENV{PETSC_ARCH} AND NOT "$ENV{PETSC_ARCH}" STREQUAL "") + list(APPEND _slepc_include_hints "${_root}/$ENV{PETSC_ARCH}/include/slepc") + list(APPEND _slepc_library_hints "${_root}/$ENV{PETSC_ARCH}/lib") + endif() +endforeach() + +if(DEFINED SLEPC_INCLUDE_DIR AND NOT "${SLEPC_INCLUDE_DIR}" MATCHES "-NOTFOUND$") + list(APPEND _slepc_include_hints ${SLEPC_INCLUDE_DIR}) +endif() +if(DEFINED SLEPC_LIBRARY AND NOT "${SLEPC_LIBRARY}" MATCHES "-NOTFOUND$") + foreach(_library IN LISTS SLEPC_LIBRARY) + if(EXISTS "${_library}") + get_filename_component(_library_dir "${_library}" DIRECTORY) + get_filename_component(_prefix_dir "${_library_dir}" DIRECTORY) + list(APPEND _slepc_library_hints "${_library_dir}") + list(APPEND _slepc_include_hints "${_prefix_dir}/include") + endif() + endforeach() +endif() +list(REMOVE_DUPLICATES _slepc_include_hints) +list(REMOVE_DUPLICATES _slepc_library_hints) + +find_package(PETSc CONFIG QUIET PATHS ${_petsc_config_paths}) +find_package(SLEPc CONFIG QUIET PATHS ${_slepc_config_paths}) +find_package(PkgConfig QUIET) +if(PkgConfig_FOUND) + pkg_check_modules(PC_PETSC QUIET IMPORTED_TARGET PETSc) + pkg_check_modules(PC_SLEPC QUIET IMPORTED_TARGET SLEPc) + if(PC_PETSC_FOUND) + list(APPEND _petsc_include_hints ${PC_PETSC_INCLUDE_DIRS}) + list(APPEND _petsc_library_hints ${PC_PETSC_LIBRARY_DIRS}) + endif() + if(PC_SLEPC_FOUND) + list(APPEND _slepc_include_hints ${PC_SLEPC_INCLUDE_DIRS}) + list(APPEND _slepc_library_hints ${PC_SLEPC_LIBRARY_DIRS}) + endif() +endif() + +function(_opensn_append_target_property output target property) + if(TARGET "${target}") + get_target_property(_property_value "${target}" "${property}") + if(_property_value AND NOT _property_value STREQUAL "_property_value-NOTFOUND") + set(${output} ${${output}} ${_property_value} PARENT_SCOPE) + endif() + endif() +endfunction() + +function(_opensn_find_imported_target output) + foreach(_target IN LISTS ARGN) + if(TARGET "${_target}") + set(${output} "${_target}" PARENT_SCOPE) + return() + endif() + endforeach() + set(${output} "" PARENT_SCOPE) +endfunction() + +_opensn_find_imported_target(_petsc_target PETSc::PETSc petsc PkgConfig::PC_PETSC) +_opensn_find_imported_target(_slepc_target SLEPc::SLEPc SLEPc::slepc slepc PkgConfig::PC_SLEPC) + +if(_petsc_target) + set(PETSC_LIBRARY "${_petsc_target}") + _opensn_append_target_property(PETSC_INCLUDE_DIR "${_petsc_target}" INTERFACE_INCLUDE_DIRECTORIES) + _opensn_append_target_property(PETSC_LIBRARY "${_petsc_target}" INTERFACE_LINK_LIBRARIES) + get_target_property(_petsc_version "${_petsc_target}" INTERFACE_VERSION) + if(_petsc_version AND NOT _petsc_version STREQUAL "_petsc_version-NOTFOUND") + set(PETSC_VERSION "${_petsc_version}") + endif() +else() + find_path(PETSC_BASE_INCLUDE_DIR + NAMES petsc.h + HINTS ${_petsc_include_hints} + DOC "PETSc base include directory" + ) + find_library(PETSC_BASE_LIBRARY + NAMES petsc + HINTS ${_petsc_library_hints} + DOC "PETSc library" + ) + set(PETSC_INCLUDE_DIR ${PETSC_BASE_INCLUDE_DIR}) + set(PETSC_LIBRARY ${PETSC_BASE_LIBRARY}) +endif() + +foreach(_include_dir IN LISTS PETSC_INCLUDE_DIR) + if(IS_DIRECTORY "${_include_dir}") + list(APPEND _petsc_include_hints "${_include_dir}") + list(APPEND _slepc_include_hints "${_include_dir}") + if(IS_DIRECTORY "${_include_dir}/slepc") + list(APPEND _slepc_include_hints "${_include_dir}/slepc") + endif() + endif() +endforeach() +list(REMOVE_DUPLICATES _petsc_include_hints) +list(REMOVE_DUPLICATES _slepc_include_hints) + +find_path(PETSC_ARCH_INCLUDE_DIR + NAMES petscconf.h + HINTS ${_petsc_include_hints} + DOC "PETSc architecture include directory" +) +find_file(PETSC_CONF_HEADER + NAMES petscconf.h + HINTS ${_petsc_include_hints} + DOC "PETSc configuration header" +) +find_file(PETSCVERSION_H + NAMES petscversion.h + HINTS ${_petsc_include_hints} + DOC "PETSc version header" ) -find_library( - PETSC_LIBRARY - petsc - PATHS - $ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib - $ENV{PETSC_DIR}/lib - $ENV{PETSC_ROOT}/$ENV{PETSC_ARCH}/lib - $ENV{PETSC_ROOT}/lib +list(APPEND PETSC_INCLUDE_DIR ${PETSC_ARCH_INCLUDE_DIR}) + +if(NOT DEFINED PETSC_VERSION) + set(PETSC_VERSION "unknown") + if(PETSCVERSION_H) + file(READ "${PETSCVERSION_H}" _petscver) + string(REGEX MATCH "#[ \t]*define[ \t]+PETSC_VERSION_MAJOR[ \t]+([0-9]+)" _major_match + "${_petscver}") + set(_petsc_major "${CMAKE_MATCH_1}") + string(REGEX MATCH "#[ \t]*define[ \t]+PETSC_VERSION_MINOR[ \t]+([0-9]+)" _minor_match + "${_petscver}") + set(_petsc_minor "${CMAKE_MATCH_1}") + string(REGEX MATCH "#[ \t]*define[ \t]+PETSC_VERSION_SUBMINOR[ \t]+([0-9]+)" _patch_match + "${_petscver}") + set(_petsc_patch "${CMAKE_MATCH_1}") + if(_major_match AND _minor_match AND _patch_match) + set(PETSC_VERSION "${_petsc_major}.${_petsc_minor}.${_petsc_patch}") + endif() + endif() +endif() + +foreach(_library IN LISTS PETSC_LIBRARY) + if(EXISTS "${_library}") + get_filename_component(_library_dir "${_library}" DIRECTORY) + get_filename_component(_prefix_dir "${_library_dir}" DIRECTORY) + list(APPEND _slepc_library_hints "${_library_dir}") + list(APPEND _slepc_include_hints "${_prefix_dir}/include") + if(IS_DIRECTORY "${_prefix_dir}/include/slepc") + list(APPEND _slepc_include_hints "${_prefix_dir}/include/slepc") + endif() + endif() +endforeach() +list(REMOVE_DUPLICATES _slepc_include_hints) +list(REMOVE_DUPLICATES _slepc_library_hints) + +if(_slepc_target) + set(SLEPC_LIBRARY "${_slepc_target}") + _opensn_append_target_property(SLEPC_INCLUDE_DIR "${_slepc_target}" INTERFACE_INCLUDE_DIRECTORIES) + _opensn_append_target_property(SLEPC_LIBRARY "${_slepc_target}" INTERFACE_LINK_LIBRARIES) +else() + find_path(SLEPC_INCLUDE_DIR + NAMES slepc.h + HINTS ${_slepc_include_hints} + DOC "SLEPc include directory" + ) + find_library(SLEPC_LIBRARY + NAMES slepc + HINTS ${_slepc_library_hints} + DOC "SLEPc library" + ) +endif() + +find_file(SLEPCVERSION_H + NAMES slepcversion.h + HINTS ${_slepc_include_hints} + DOC "SLEPc version header" ) -set(PETSC_VERSION "unknown") -find_file(PETSCVERSION_H petscversion.h - PATHS - $ENV{PETSC_DIR}/include - $ENV{PETSC_ROOT}/include +set(SLEPC_VERSION "unknown") +if(SLEPCVERSION_H) + file(READ "${SLEPCVERSION_H}" _slepcver) + string(REGEX MATCH "#[ \t]*define[ \t]+SLEPC_VERSION_MAJOR[ \t]+([0-9]+)" _major_match + "${_slepcver}") + set(_slepc_major "${CMAKE_MATCH_1}") + string(REGEX MATCH "#[ \t]*define[ \t]+SLEPC_VERSION_MINOR[ \t]+([0-9]+)" _minor_match + "${_slepcver}") + set(_slepc_minor "${CMAKE_MATCH_1}") + string(REGEX MATCH "#[ \t]*define[ \t]+SLEPC_VERSION_SUBMINOR[ \t]+([0-9]+)" _patch_match + "${_slepcver}") + set(_slepc_patch "${CMAKE_MATCH_1}") + if(_major_match AND _minor_match AND _patch_match) + set(SLEPC_VERSION "${_slepc_major}.${_slepc_minor}.${_slepc_patch}") + endif() +endif() + +list(APPEND PETSC_INCLUDE_DIR ${SLEPC_INCLUDE_DIR}) +list(APPEND PETSC_LIBRARY ${SLEPC_LIBRARY}) +list(REMOVE_DUPLICATES PETSC_INCLUDE_DIR) +list(REMOVE_DUPLICATES PETSC_LIBRARY) + +find_package_handle_standard_args( + SLEPc + REQUIRED_VARS SLEPC_LIBRARY SLEPC_INCLUDE_DIR + VERSION_VAR SLEPC_VERSION + NAME_MISMATCHED ) -mark_as_advanced(PETSCVERSION_H) -if (PETSCVERSION_H) - file(READ ${PETSCVERSION_H} PETSC_VERSION_FILE) - string(REGEX MATCH "define[ ]+PETSC_VERSION_MAJOR[ ]+([0-9]+)" TMP "${PETSC_VERSION_FILE}") - set(PETSC_VERSION_MAJOR ${CMAKE_MATCH_1}) - string(REGEX MATCH "define[ ]+PETSC_VERSION_MINOR[ ]+([0-9]+)" TMP "${PETSC_VERSION_FILE}") - set(PETSC_VERSION_MINOR ${CMAKE_MATCH_1}) - string(REGEX MATCH "define[ ]+PETSC_VERSION_SUBMINOR[ ]+([0-9]+)" TMP "${PETSC_VERSION_FILE}") - set(PETSC_VERSION_PATCH ${CMAKE_MATCH_1}) - set(PETSC_VERSION "${PETSC_VERSION_MAJOR}.${PETSC_VERSION_MINOR}.${PETSC_VERSION_PATCH}") +set(SLEPC_FOUND ${SLEPc_FOUND}) + +if(NOT SLEPc_FOUND) + # OpenSn requires SLEPc wherever PETSc is required. + set(PETSC_INCLUDE_DIR PETSC_INCLUDE_DIR-NOTFOUND) + set(PETSC_LIBRARY PETSC_LIBRARY-NOTFOUND) + if(PETSc_FIND_REQUIRED) + message(FATAL_ERROR "PETSc/SLEPc support is required, but SLEPc was not found.") + endif() endif() -include(FindPackageHandleStandardArgs) find_package_handle_standard_args( - PETSc - REQUIRED_VARS PETSC_LIBRARY PETSC_INCLUDE_DIR - VERSION_VAR PETSC_VERSION + PETSc + REQUIRED_VARS PETSC_LIBRARY PETSC_INCLUDE_DIR PETSC_CONF_HEADER + VERSION_VAR PETSC_VERSION ) +set(PETSC_FOUND ${PETSc_FOUND}) diff --git a/framework/math/linear_solver/linear_eigen_solver.h b/framework/math/linear_solver/linear_eigen_solver.h index d23c8a3bf6..1d40acc200 100644 --- a/framework/math/linear_solver/linear_eigen_solver.h +++ b/framework/math/linear_solver/linear_eigen_solver.h @@ -14,7 +14,10 @@ class LinearEigenSolver : public LinearSolver public: enum class IterativeMethod : int { - POWER + KRYLOV_SCHUR, + POWER, + ARNOLDI, + SUBSPACE }; LinearEigenSolver(IterativeMethod method, std::shared_ptr ctx) @@ -26,8 +29,14 @@ class LinearEigenSolver : public LinearSolver { switch (method_) { + case IterativeMethod::KRYLOV_SCHUR: + return "KRYLOV-SCHUR"; case IterativeMethod::POWER: return "POWER ITERATION"; + case IterativeMethod::ARNOLDI: + return "ARNOLDI"; + case IterativeMethod::SUBSPACE: + return "SUBSPACE"; default: return "UNKNOWN ITERATIVE METHOD"; } diff --git a/framework/math/linear_solver/slepc_linear_eigen_solver.cc b/framework/math/linear_solver/slepc_linear_eigen_solver.cc new file mode 100644 index 0000000000..b85ed5d358 --- /dev/null +++ b/framework/math/linear_solver/slepc_linear_eigen_solver.cc @@ -0,0 +1,103 @@ +// SPDX-FileCopyrightText: 2025 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#include "framework/math/linear_solver/slepc_linear_eigen_solver.h" +#include "framework/runtime.h" + +namespace opensn +{ + +SLEPcLinearEigenSolver::SLEPcLinearEigenSolver(IterativeMethod method, + std::shared_ptr context_ptr) + : LinearEigenSolver(method, context_ptr), + A_(nullptr), + x_(nullptr), + eps_(nullptr), + eps_type_("krylovschur"), + num_local_dofs_(0), + num_global_dofs_(0), + system_set_(false) +{ +} + +SLEPcLinearEigenSolver::~SLEPcLinearEigenSolver() +{ + MatDestroy(&A_); + VecDestroy(&x_); + EPSDestroy(&eps_); +} + +void +SLEPcLinearEigenSolver::PreSetupCallback() +{ +} + +void +SLEPcLinearEigenSolver::SetOptions() +{ +} + +void +SLEPcLinearEigenSolver::SetSolverContext() +{ +} + +void +SLEPcLinearEigenSolver::SetConvergenceTest() +{ +} + +void +SLEPcLinearEigenSolver::SetMonitor() +{ +} + +void +SLEPcLinearEigenSolver::SetPreconditioner() +{ +} + +void +SLEPcLinearEigenSolver::PostSetupCallback() +{ +} + +void +SLEPcLinearEigenSolver::Setup() +{ + if (IsSystemSet()) + return; + + PreSetupCallback(); + EPSCreate(opensn::mpi_comm, &eps_); + SetOptions(); + SetSolverContext(); + SetConvergenceTest(); + SetMonitor(); + SetSystemSize(); + SetSystem(); + SetPreconditioner(); + PostSetupCallback(); + + system_set_ = true; +} + +void +SLEPcLinearEigenSolver::PreSolveCallback() +{ +} + +void +SLEPcLinearEigenSolver::PostSolveCallback() +{ +} + +void +SLEPcLinearEigenSolver::Solve() +{ + PreSolveCallback(); + SetInitialGuess(); + PostSolveCallback(); +} + +} // namespace opensn diff --git a/framework/math/linear_solver/slepc_linear_eigen_solver.h b/framework/math/linear_solver/slepc_linear_eigen_solver.h new file mode 100644 index 0000000000..963489010c --- /dev/null +++ b/framework/math/linear_solver/slepc_linear_eigen_solver.h @@ -0,0 +1,82 @@ +// SPDX-FileCopyrightText: 2024 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#pragma once + +#include "framework/math/linear_solver/linear_eigen_solver.h" +#include "framework/math/linear_solver/linear_solver_context.h" +#include +#include +#include +#include + +namespace opensn +{ + +class SLEPcLinearEigenSolver : public LinearEigenSolver +{ +public: + struct ToleranceOptions + { + double residual_absolute = 1.0e-6; + int maximum_iterations = 100; + } tolerance_options; + + SLEPcLinearEigenSolver(IterativeMethod method, std::shared_ptr context_ptr); + + ~SLEPcLinearEigenSolver() override; + + ToleranceOptions& GetToleranceOptions() { return tolerance_options; } + + std::string GetEPSType() { return eps_type_; } + + void SetEPSType(const std::string& type) + { + eps_type_ = type; + if (type == "power") + method_ = IterativeMethod::POWER; + else if (type == "subspace") + method_ = IterativeMethod::SUBSPACE; + else if (type == "arnoldi") + method_ = IterativeMethod::ARNOLDI; + else if (type == "krylovschur") + method_ = IterativeMethod::KRYLOV_SCHUR; + else + throw std::runtime_error("Invalid iterative method type for linear eigen solver"); + } + + void ApplyToleranceOptions(); + + /// Set up the linaer solver + void Setup() override; + + /// Solve the system + void Solve() override; + +protected: + bool IsSystemSet() const { return system_set_; } + virtual void PreSetupCallback(); + virtual void SetOptions(); + virtual void SetSolverContext(); + virtual void SetConvergenceTest(); + virtual void SetMonitor(); + virtual void SetPreconditioner(); + virtual void SetSystemSize() = 0; + virtual void SetSystem() = 0; + virtual void PostSetupCallback(); + virtual void PreSolveCallback(); + virtual void SetInitialGuess() = 0; + virtual void PostSolveCallback(); + + Mat A_; + Vec x_; + EPS eps_; + std::string eps_type_; + int64_t num_local_dofs_; + int64_t num_global_dofs_; + +private: + bool system_set_; +}; + +} // namespace opensn diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/slepc_linear_keigen_solver.cc b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/slepc_linear_keigen_solver.cc new file mode 100644 index 0000000000..f08952c6c1 --- /dev/null +++ b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/slepc_linear_keigen_solver.cc @@ -0,0 +1,373 @@ +// SPDX-FileCopyrightText: 2025 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/slepc_linear_keigen_solver.h" +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/ags_linear_solver.h" +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/wgs_context.h" +#include "modules/linear_boltzmann_solvers/lbs_problem/compute/lbs_compute.h" +#include "modules/linear_boltzmann_solvers/lbs_problem/vecops/lbs_vecops.h" +#include "framework/math/petsc_utils/petsc_utils.h" +#include "framework/math/math.h" +#include "framework/logging/log.h" +#include "framework/utils/timer.h" +#include "framework/runtime.h" +#include +#include +#include +#include +#include +#include +#include + +namespace opensn +{ + +namespace +{ + +struct SavedWGSSourceScopes +{ + std::vector contexts; + std::vector lhs_scopes; + std::vector rhs_scopes; + + explicit SavedWGSSourceScopes(DiscreteOrdinatesProblem& do_problem) + { + const auto num_wgs = do_problem.GetNumWGSSolvers(); + contexts.reserve(num_wgs); + lhs_scopes.reserve(num_wgs); + rhs_scopes.reserve(num_wgs); + + for (size_t i = 0; i < num_wgs; ++i) + { + auto wgs_solver = do_problem.GetWGSSolver(i); + OpenSnLogicalErrorIf(not wgs_solver, + do_problem.GetName() + ": Null WGS solver in SLEPc k-eigen solve."); + + auto context = wgs_solver->GetContext(); + auto wgs_context = std::dynamic_pointer_cast(context); + OpenSnLogicalErrorIf(not wgs_context, do_problem.GetName() + ": Cast to WGSContext failed."); + + contexts.push_back(wgs_context.get()); + lhs_scopes.push_back(wgs_context->lhs_src_scope); + rhs_scopes.push_back(wgs_context->rhs_src_scope); + + wgs_context->lhs_src_scope.Unset(APPLY_WGS_FISSION_SOURCES); + wgs_context->rhs_src_scope.Unset(APPLY_AGS_FISSION_SOURCES); + } + } + + ~SavedWGSSourceScopes() + { + for (size_t i = 0; i < contexts.size(); ++i) + { + contexts[i]->lhs_src_scope = lhs_scopes[i]; + contexts[i]->rhs_src_scope = rhs_scopes[i]; + } + } +}; + +struct SavedTransportState +{ + struct AngleAggregationState + { + AngleAggregation* angle_agg = nullptr; + std::vector psi_old; + std::vector psi_new; + }; + + DiscreteOrdinatesProblem& do_problem; + std::vector phi_old; + std::vector phi_new; + std::vector q_moments; + std::vector angle_agg_states; + + explicit SavedTransportState(DiscreteOrdinatesProblem& do_problem) + : do_problem(do_problem), + phi_old(do_problem.GetPhiOldLocal()), + phi_new(do_problem.GetPhiNewLocal()), + q_moments(do_problem.GetQMomentsLocal()) + { + angle_agg_states.reserve(do_problem.GetGroupsets().size()); + for (const auto& groupset : do_problem.GetGroupsets()) + { + if (not groupset.angle_agg) + continue; + angle_agg_states.push_back({groupset.angle_agg.get(), + groupset.angle_agg->GetOldDelayedAngularDOFsAsSTLVector(), + groupset.angle_agg->GetNewDelayedAngularDOFsAsSTLVector()}); + } + } + + ~SavedTransportState() noexcept + { + try + { + Restore(); + } + catch (...) + { + std::terminate(); + } + } + + void Restore() + { + CopyInto(do_problem.GetPhiOldLocal(), phi_old); + CopyInto(do_problem.GetPhiNewLocal(), phi_new); + CopyInto(do_problem.GetQMomentsLocal(), q_moments); + for (const auto& state : angle_agg_states) + { + state.angle_agg->SetOldDelayedAngularDOFsFromSTLVector(state.psi_old); + state.angle_agg->SetNewDelayedAngularDOFsFromSTLVector(state.psi_new); + } + } + + static void CopyInto(std::vector& destination, const std::vector& source) + { + assert(destination.size() == source.size() && "SavedTransportState restore size mismatch."); + std::copy(source.begin(), source.end(), destination.begin()); + } +}; + +PetscErrorCode +ApplyLinearKEigenOperator(SLEPcLinearKEigenContext& ctx, Vec x, Vec y) +{ + auto& do_problem = ctx.do_problem; + auto& phi_old_local = do_problem->GetPhiOldLocal(); + auto& phi_new_local = do_problem->GetPhiNewLocal(); + auto& q_moments_local = do_problem->GetQMomentsLocal(); + auto active_set_source_function = do_problem->GetActiveSetSourceFunction(); + auto ags_solver = do_problem->GetAGSSolver(); + + LBSVecOps::SetPrimarySTLvectorFromGroupScopedPETScVec( + *do_problem, 0, do_problem->GetNumGroups() - 1, x, PhiSTLOption::PHI_OLD); + + Set(q_moments_local, 0.0); + + // Apply the linear k-eigenvalue operator: + // A phi = (I - D L^-1 M S)^-1 D L^-1 F phi. + // The fission source is the explicit input to the operator. The AGS solve applies the + // scattering inverse with fission removed from the WGS source scopes. + for (const auto& groupset : do_problem->GetGroupsets()) + active_set_source_function(groupset, + q_moments_local, + phi_old_local, + APPLY_AGS_FISSION_SOURCES | APPLY_WGS_FISSION_SOURCES); + + OpenSnLogicalErrorIf(not ags_solver, "SLEPcLinearKEigenSolver: AGS solver is not available."); + ags_solver->Solve(); + + LBSVecOps::SetGroupScopedPETScVecFromPrimarySTLvector( + *do_problem, 0, do_problem->GetNumGroups() - 1, y, phi_new_local); + + return 0; +} + +// Shell matrix multiplication: y = A * x +PetscErrorCode +ShellMult(Mat M, Vec x, Vec y) +{ + void* raw_ctx = nullptr; + MatShellGetContext(M, static_cast(&raw_ctx)); + + auto* ctx = static_cast(raw_ctx); + SavedTransportState saved_transport_state(*ctx->do_problem); + + return ApplyLinearKEigenOperator(*ctx, x, y); +} + +} // namespace + +static PetscErrorCode +SLEPcLinearKEigenMonitor(EPS unused_eps, + PetscInt its, + PetscInt unused_nconv, + PetscScalar* unused_eigr, // NOLINT(readability-non-const-parameter) + PetscScalar* unused_eigi, // NOLINT(readability-non-const-parameter) + PetscReal* errest, + PetscInt nest, + void* ctx) +{ + (void)unused_eps; + (void)unused_nconv; + (void)unused_eigr; + (void)unused_eigi; + + const bool has_error_estimate = nest > 0 and errest != nullptr; + if (ctx) + { + auto* kctx = static_cast(ctx); + if (has_error_estimate) + kctx->last_eps_residual = errest[0]; + } + std::stringstream iter_info; + iter_info << program_timer.GetTimeString() << " SLEPc EPS iteration = " << its; + if (has_error_estimate) + iter_info << ", residual = " << std::scientific << std::setprecision(6) << errest[0] + << std::defaultfloat; + else + iter_info << ", status = iterating"; + log.Log() << iter_info.str(); + + return 0; +} + +void +SLEPcLinearKEigenSolver::SetMonitor() +{ +} + +void +SLEPcLinearKEigenSolver::SetSystemSize() +{ + auto ctx = std::dynamic_pointer_cast(context_ptr_); + // SLEPc solves the fission-production eigenproblem in scalar-flux space. Delayed angular + // fluxes are still solved as part of each transport inverse application, but they are internal + // operator state rather than EPS eigenvector entries. + auto sizes = ctx->do_problem->LBSProblem::GetNumPhiIterativeUnknowns(); + num_local_dofs_ = static_cast(sizes.first); + num_global_dofs_ = static_cast(sizes.second); +} + +void +SLEPcLinearKEigenSolver::SetSystem() +{ + // Create shell matrix A + MatCreateShell(PETSC_COMM_WORLD, + static_cast(num_local_dofs_), + static_cast(num_local_dofs_), + static_cast(num_global_dofs_), + static_cast(num_global_dofs_), + context_ptr_.get(), + &A_); + // PETSc erases the shell-operation type to PetscErrorCodeFn*, so a cast is required here. + // NOLINTNEXTLINE(cppcoreguidelines-pro-type-cstyle-cast,modernize-redundant-void-arg) + MatShellSetOperation(A_, MATOP_MULT, (void (*)(void))ShellMult); + + // Create x + x_ = CreateVector(num_local_dofs_, num_global_dofs_); +} + +void +SLEPcLinearKEigenSolver::SetInitialGuess() +{ + auto ctx = std::dynamic_pointer_cast(context_ptr_); + LBSVecOps::SetGroupScopedPETScVecFromPrimarySTLvector(*ctx->do_problem, + 0, + ctx->do_problem->GetNumGroups() - 1, + x_, + ctx->do_problem->GetPhiOldLocal()); +} + +void +SLEPcLinearKEigenSolver::Solve() +{ + PreSolveCallback(); + auto ctx = std::dynamic_pointer_cast(context_ptr_); + SavedWGSSourceScopes saved_wgs_scopes(*ctx->do_problem); + SetInitialGuess(); + + EPSSetFromOptions(eps_); + EPSSetOperators(eps_, A_, nullptr); + EPSSetProblemType(eps_, EPS_NHEP); + EPSSetWhichEigenpairs(eps_, EPS_LARGEST_MAGNITUDE); + EPSSetTolerances(eps_, tolerance_options.residual_absolute, tolerance_options.maximum_iterations); + EPSSetType(eps_, eps_type_.c_str()); + EPSSetInitialSpace(eps_, 1, &x_); + EPSMonitorSet(eps_, SLEPcLinearKEigenMonitor, context_ptr_.get(), nullptr); + EPSSolve(eps_); + + // Check for convergence + PetscInt nconv = 0; + EPSGetConverged(eps_, &nconv); + if (nconv == 0) + { + OpenSnLogicalError(program_timer.GetTimeString() + " SLEPc EPS: No eigenpairs converged."); + } + + PostSolveCallback(); +} + +void +SLEPcLinearKEigenSolver::PreSolveCallback() +{ + log.Log() << "Executing SLEPc Eigenvalue Problem Solver with iterative method " + << GetIterativeMethodName() << std::endl; + + auto ctx = std::dynamic_pointer_cast(context_ptr_); + ctx->last_eps_residual = 1.0; + ctx->last_operator_residual = 1.0; +} + +void +SLEPcLinearKEigenSolver::PostSolveCallback() +{ + auto ctx = std::dynamic_pointer_cast(context_ptr_); + auto& do_problem = ctx->do_problem; + auto& phi_old_local = do_problem->GetPhiOldLocal(); + auto& phi_new_local = do_problem->GetPhiNewLocal(); + + // Extract the converged eigenvector and compute k as the fission-production ratio. + PetscScalar kr = 0.0; + PetscScalar ki = 0.0; + Vec x = CreateVector(num_local_dofs_, num_global_dofs_); + EPSGetEigenpair(eps_, 0, &kr, &ki, x, nullptr); + LBSVecOps::SetPrimarySTLvectorFromGroupScopedPETScVec( + *do_problem, 0, do_problem->GetNumGroups() - 1, x, PhiSTLOption::PHI_OLD); + double F_prev = ComputeFissionProduction(*do_problem, phi_old_local); + + // Perform a transport solve, get the resulting eigenvector as phi_new, and compute fission rate + Vec y = CreateVector(num_local_dofs_, num_global_dofs_); + ApplyLinearKEigenOperator(*ctx, x, y); + LBSVecOps::SetPrimarySTLvectorFromGroupScopedPETScVec( + *do_problem, 0, do_problem->GetNumGroups() - 1, y, PhiSTLOption::PHI_NEW); + double F = ComputeFissionProduction(*do_problem, phi_new_local); + + // Compute k-eigenvalue + OpenSnLogicalErrorIf(F_prev == 0.0, + "SLEPcLinearKEigenSolver: previous fission production is zero. " + "Cannot compute k-eigenvalue."); + ctx->eigenvalue = F / F_prev; + log.Log() << "Final k-eigenvalue: " << std::setprecision(7) << ctx->eigenvalue << "\n"; + + // Compute operator residual ||A x - k x|| / ||x|| using the same operator SLEPc solved. + Vec residual = CreateVector(num_local_dofs_, num_global_dofs_); + OpenSnPETScCall(VecCopy(y, residual)); + OpenSnPETScCall(VecAXPY(residual, -ctx->eigenvalue, x)); + PetscReal x_norm = 0.0; + PetscReal residual_norm = 0.0; + OpenSnPETScCall(VecNorm(x, NORM_2, &x_norm)); + OpenSnPETScCall(VecNorm(residual, NORM_2, &residual_norm)); + ctx->last_operator_residual = (x_norm > 0.0) ? residual_norm / x_norm : residual_norm; + log.Log() << "Operator residual: " << std::setprecision(7) << ctx->last_operator_residual << "\n"; + + OpenSnLogicalErrorIf(not std::isfinite(ctx->last_operator_residual), + "SLEPcLinearKEigenSolver: final operator residual is not finite."); + + const auto ags_solver = do_problem->GetAGSSolver(); + const double residual_tolerance = std::max({10.0 * tolerance_options.residual_absolute, + ags_solver ? 10.0 * ags_solver->GetTolerance() : 0.0, + 1.0e-10}); + if (ctx->last_operator_residual > residual_tolerance) + { + std::stringstream warning; + warning << "SLEPcLinearKEigenSolver: final operator residual " << std::scientific + << std::setprecision(6) << ctx->last_operator_residual + << " exceeds diagnostic tolerance " << residual_tolerance + << ". This residual includes the accuracy of the inner transport solve."; + log.Log0Warning() << warning.str(); + } + + // Leave the problem state consistent with the operator application. + LBSVecOps::SetPrimarySTLvectorFromGroupScopedPETScVec( + *do_problem, 0, do_problem->GetNumGroups() - 1, x, PhiSTLOption::PHI_OLD); + LBSVecOps::SetPrimarySTLvectorFromGroupScopedPETScVec( + *do_problem, 0, do_problem->GetNumGroups() - 1, y, PhiSTLOption::PHI_NEW); + + VecDestroy(&x); + VecDestroy(&y); + VecDestroy(&residual); +} + +} // namespace opensn diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/slepc_linear_keigen_solver.h b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/slepc_linear_keigen_solver.h new file mode 100644 index 0000000000..91306002ad --- /dev/null +++ b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/slepc_linear_keigen_solver.h @@ -0,0 +1,52 @@ +// SPDX-FileCopyrightText: 2025 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#pragma once + +#include "framework/math/linear_solver/slepc_linear_eigen_solver.h" +#include "framework/math/linear_solver/linear_solver_context.h" +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.h" +#include + +namespace opensn +{ + +/// Context for the linear k-eigenvalue solve +struct SLEPcLinearKEigenContext : public LinearEigenContext +{ + std::shared_ptr do_problem; + double last_eps_residual = 1.0; + double last_operator_residual = 1.0; + + explicit SLEPcLinearKEigenContext(std::shared_ptr do_problem) + : do_problem(std::move(do_problem)) + { + eigenvalue = 1.0; + } + + ~SLEPcLinearKEigenContext() override = default; +}; + +/// A SLEPc shell-based linear k-eigenvalue solver +class SLEPcLinearKEigenSolver : public SLEPcLinearEigenSolver +{ +public: + explicit SLEPcLinearKEigenSolver(std::shared_ptr context) + : SLEPcLinearEigenSolver(IterativeMethod::KRYLOV_SCHUR, std::move(context)) + { + } + + ~SLEPcLinearKEigenSolver() override = default; + + void Solve() override; + +protected: + void SetMonitor() override; + void SetSystemSize() override; + void SetSystem() override; + void SetInitialGuess() override; + void PostSolveCallback() override; + void PreSolveCallback() override; +}; + +} // namespace opensn diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/solvers/slepc_keigen_solver.cc b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/solvers/slepc_keigen_solver.cc new file mode 100644 index 0000000000..84c22cb6bd --- /dev/null +++ b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/solvers/slepc_keigen_solver.cc @@ -0,0 +1,89 @@ +// SPDX-FileCopyrightText: 2025 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/solvers/slepc_keigen_solver.h" +#include "modules/linear_boltzmann_solvers/lbs_problem/compute/lbs_compute.h" +#include "modules/linear_boltzmann_solvers/lbs_problem/vecops/lbs_vecops.h" +#include "framework/object_factory.h" +#include "framework/logging/log.h" +#include "framework/runtime.h" + +namespace opensn +{ + +OpenSnRegisterObjectInNamespace(lbs, SLEPcKEigenSolver); + +InputParameters +SLEPcKEigenSolver::GetInputParameters() +{ + InputParameters params = Solver::GetInputParameters(); + + params.SetGeneralDescription("Implementation of a linear k-eigenvalue solver using SLEPc EPS"); + params.ChangeExistingParamToOptional("name", "SLEPcEigenValueProblemSolver"); + + params.AddRequiredParameter>("lbs_problem", + "An existing discrete ordinates problem"); + params.AddOptionalParameter("max_iters", 50, "Maximum iterations"); + params.AddOptionalParameter("eigen_tol", 1.0e-8, "Absolute tolerance"); + params.AddOptionalParameter("reset_phi0", true, "If true, reinitializes scalar fluxes to 1.0"); + params.AddOptionalParameter("eps_type", "krylovschur", "EPS solver type"); + params.ConstrainParameterRange( + "eps_type", AllowableRangeList::New({"power", "subspace", "arnoldi", "krylovschur"})); + + return params; +} + +std::shared_ptr +SLEPcKEigenSolver::Create(const ParameterBlock& params) +{ + auto& factory = opensn::ObjectFactory::GetInstance(); + return factory.Create("lbs::SLEPcKEigenSolver", params); +} + +SLEPcKEigenSolver::SLEPcKEigenSolver(const InputParameters& params) + : Solver(params), + do_problem_(params.GetSharedPtrParam("lbs_problem")), + context_(std::make_shared(do_problem_)), + solver_(context_), + reset_phi0_(params.GetParamValue("reset_phi0")) +{ + auto& tolerances = solver_.GetToleranceOptions(); + tolerances.residual_absolute = params.GetParamValue("eigen_tol"); + tolerances.maximum_iterations = params.GetParamValue("max_iters"); + auto type = params.GetParamValue("eps_type"); + solver_.SetEPSType(type); +} + +void +SLEPcKEigenSolver::Initialize() +{ + OpenSnInvalidArgumentIf(do_problem_->IsTimeDependent(), + GetName() + ": Problem is in time-dependent mode. Call problem." + "SetSteadyStateMode() before initializing this solver."); + + if (reset_phi0_) + LBSVecOps::SetPhiVectorScalarValues(*do_problem_, PhiSTLOption::PHI_OLD, 1.0); +} + +void +SLEPcKEigenSolver::Execute() +{ + solver_.Setup(); + solver_.Solve(); + + if (do_problem_->GetOptions().use_precursors) + { + ComputePrecursors(*do_problem_); + do_problem_->ScalePrecursors(1.0 / context_->eigenvalue); + } + + log.Log() << "SLEPcKEigenSolver execution completed\n\n"; +} + +double +SLEPcKEigenSolver::GetEigenvalue() const +{ + return context_->eigenvalue; +} + +} // namespace opensn diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/solvers/slepc_keigen_solver.h b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/solvers/slepc_keigen_solver.h new file mode 100644 index 0000000000..47ad246259 --- /dev/null +++ b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/solvers/slepc_keigen_solver.h @@ -0,0 +1,36 @@ +// SPDX-FileCopyrightText: 2025 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#pragma once + +#include "modules/solver.h" +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/slepc_linear_keigen_solver.h" +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.h" + +namespace opensn +{ + +class SLEPcKEigenSolver : public Solver +{ +private: + std::shared_ptr do_problem_; + std::shared_ptr context_; + SLEPcLinearKEigenSolver solver_; + + bool reset_phi0_; + +public: + explicit SLEPcKEigenSolver(const InputParameters& params); + + void Initialize() override; + void Execute() override; + + /// Return the k-eigenvalue + double GetEigenvalue() const; + +public: + static InputParameters GetInputParameters(); + static std::shared_ptr Create(const ParameterBlock& params); +}; + +} // namespace opensn diff --git a/python/lib/py_app.cc b/python/lib/py_app.cc index a950009cef..1bf1b57e22 100644 --- a/python/lib/py_app.cc +++ b/python/lib/py_app.cc @@ -63,6 +63,7 @@ PyApp::PyApp(const mpi::Communicator& comm) Console::BindModule(WrapLBS); Console::BindModule(WrapSteadyState); Console::BindModule(WrapTransient); + Console::BindModule(WrapSLEPcKEigen); Console::BindModule(WrapNLKEigen); Console::BindModule(WrapPIteration); Console::BindModule(WrapDiscreteOrdinatesKEigenAcceleration); diff --git a/python/lib/py_wrappers.h b/python/lib/py_wrappers.h index 32d76f13a1..3317806f84 100644 --- a/python/lib/py_wrappers.h +++ b/python/lib/py_wrappers.h @@ -227,6 +227,7 @@ void WrapSolver(py::module& slv); void WrapLBS(py::module& slv); void WrapSteadyState(py::module& slv); void WrapTransient(py::module& slv); +void WrapSLEPcKEigen(py::module& slv); void WrapNLKEigen(py::module& slv); void WrapPIteration(py::module& slv); void WrapDiscreteOrdinatesKEigenAcceleration(py::module& slv); diff --git a/python/lib/solver.cc b/python/lib/solver.cc index 6c7fc3ced5..95c487764f 100644 --- a/python/lib/solver.cc +++ b/python/lib/solver.cc @@ -15,6 +15,7 @@ #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/io/discrete_ordinates_problem_io.h" #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/solvers/transient_solver.h" #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/solvers/steady_state_solver.h" +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/solvers/slepc_keigen_solver.h" #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/solvers/nl_keigen_solver.h" #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/solvers/pi_keigen_solver.h" #include "modules/linear_boltzmann_solvers/lbs_problem/io/lbs_problem_io.h" @@ -1630,6 +1631,48 @@ WrapTransient(py::module& slv) // clang-format on } +// Wrap linear k-eigen solver +void +WrapSLEPcKEigen(py::module& slv) +{ + // clang-format off + // Linear k-eigen solver + auto linear_k_eigen_solver = py::class_, + Solver>( + slv, + "SLEPcKEigenSolver", + R"( + Linear k-eigenvalue solver. + + Wrapper of :cpp:class:`opensn::SLEPcKEigenSolver`. + )" + ); + linear_k_eigen_solver.def( + py::init( + [](py::kwargs& params) + { + return SLEPcKEigenSolver::Create(kwargs_to_param_block(params)); + } + ), + R"( + Construct a linear k-eigenvalue solver. + + Parameters + ---------- + lbs_problem: pyopensn.solver.DiscreteOrdinatesProblem + Existing discrete ordinates problem instance. + )" + ); + linear_k_eigen_solver.def( + "GetEigenvalue", + &SLEPcKEigenSolver::GetEigenvalue, + R"( + Return the current k‑eigenvalue. + )" + ); + // clang-format on +} + // Wrap non-linear k-eigen solver void WrapNLKEigen(py::module& slv) @@ -2164,6 +2207,7 @@ py_solver(py::module& pyopensn) WrapLBS(slv); WrapSteadyState(slv); WrapTransient(slv); + WrapSLEPcKEigen(slv); WrapNLKEigen(slv); WrapDiscreteOrdinatesKEigenAcceleration(slv); WrapPIteration(slv); diff --git a/python/main.cc b/python/main.cc index 67e8ded475..f11c29818a 100644 --- a/python/main.cc +++ b/python/main.cc @@ -3,18 +3,68 @@ #include "python/lib/py_app.h" #include "mpicpp-lite/mpicpp-lite.h" -#include "petsc.h" +#include "slepc.h" #include #include +#include +#include namespace mpi = mpicpp_lite; +namespace +{ + +std::vector +BuildSLEPcArgv(int argc, char** argv) +{ + std::vector slepc_argv; + slepc_argv.reserve(static_cast(argc)); + slepc_argv.push_back(argv[0]); + + for (int k = 1; k < argc; ++k) + { + const std::string arg = argv[k]; + + // These options are consumed by the OpenSn Python app, not PETSc/SLEPc. + if (arg == "-h" or arg == "--help" or arg == "-c" or arg == "--suppress-color") + continue; + + if (arg == "-v" or arg == "--verbose" or arg == "-i" or arg == "--input" or arg == "-p" or + arg == "--py") + { + ++k; // Skip the accompanying value. + continue; + } + + if (arg.starts_with("--verbose=") or arg.starts_with("--input=") or arg.starts_with("--py=") or + arg.starts_with("--caliper=")) + continue; + + if (arg == "--caliper") + { + // cxxopts uses an implicit value when no explicit value is provided. + if (k + 1 < argc and argv[k + 1][0] != '-') + ++k; + continue; + } + + slepc_argv.push_back(argv[k]); + } + + return slepc_argv; +} + +} // namespace + int main(int argc, char** argv) // NOLINT(bugprone-exception-escape) { mpi::Environment env(argc, argv, mpi::ThreadSupport::MULTIPLE); - PetscCall(PetscInitializeNoArguments()); // NOLINT(bugprone-casting-through-void) + auto slepc_argv = BuildSLEPcArgv(argc, argv); + int slepc_argc = static_cast(slepc_argv.size()); + char** slepc_argv_ptr = slepc_argv.data(); + PetscCall(SlepcInitialize(&slepc_argc, &slepc_argv_ptr, nullptr, nullptr)); int retval = EXIT_SUCCESS; try @@ -29,7 +79,7 @@ main(int argc, char** argv) // NOLINT(bugprone-exception-escape) retval = EXIT_FAILURE; } - PetscFinalize(); + SlepcFinalize(); return retval; } diff --git a/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_5g_upscatter_slepc.py b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_5g_upscatter_slepc.py new file mode 100644 index 0000000000..cf07e6a249 --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_5g_upscatter_slepc.py @@ -0,0 +1,236 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +3D 5G k-eigenvalue test with five groupsets and upscatter. + +The homogeneous fully reflecting problem has a spatially constant solution. +The reference k is the dominant eigenvalue of inv(T - S) F for the XS data in +simple_5g_upscatter_fissile.xs. +""" + +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.xs import MultiGroupXS + from pyopensn.aquad import GLCProductQuadrature3DXYZ + from pyopensn.solver import ( + DiscreteOrdinatesProblem, + PowerIterationKEigenSolver, + SLEPcKEigenSolver, + ) + + +SIGMA_T = [1.0, 1.0, 1.0, 1.0, 1.0] +SCATTER = [ + [0.20, 0.02, 0.03, 0.01, 0.02], + [0.04, 0.25, 0.02, 0.03, 0.01], + [0.06, 0.05, 0.22, 0.04, 0.02], + [0.02, 0.06, 0.07, 0.24, 0.03], + [0.01, 0.03, 0.04, 0.05, 0.21], +] +PRODUCTION = [ + [0.275, 0.330, 0.110, 0.055, 0.044], + [0.125, 0.150, 0.050, 0.025, 0.020], + [0.060, 0.072, 0.024, 0.012, 0.010], + [0.030, 0.036, 0.012, 0.006, 0.005], + [0.010, 0.012, 0.004, 0.002, 0.002], +] +EXPECTED_K = 0.6553715484430904 + + +def solve_dense(matrix, rhs): + """Small dense Gaussian solve with partial pivoting.""" + matrix = [row[:] for row in matrix] + rhs = rhs[:] + n = len(rhs) + for i in range(n): + pivot = max(range(i, n), key=lambda r: abs(matrix[r][i])) + matrix[i], matrix[pivot] = matrix[pivot], matrix[i] + rhs[i], rhs[pivot] = rhs[pivot], rhs[i] + scale = matrix[i][i] + if abs(scale) < 1.0e-14: + raise RuntimeError("Singular reference matrix.") + for j in range(i, n): + matrix[i][j] /= scale + rhs[i] /= scale + for r in range(n): + if r == i: + continue + factor = matrix[r][i] + for j in range(i, n): + matrix[r][j] -= factor * matrix[i][j] + rhs[r] -= factor * rhs[i] + return rhs + + +def apply_reference_operator(phi): + loss = [ + [(SIGMA_T[g] if g == gp else 0.0) - SCATTER[g][gp] for gp in range(5)] + for g in range(5) + ] + fission_source = [ + sum(PRODUCTION[g][gp] * phi[gp] for gp in range(5)) + for g in range(5) + ] + return solve_dense(loss, fission_source) + + +def compute_reference_k(): + phi = [1.0, 1.0, 1.0, 1.0, 1.0] + k = 1.0 + for _ in range(200): + new_phi = apply_reference_operator(phi) + k = sum(new_phi) / sum(phi) + norm = sum(new_phi) + phi = [x / norm for x in new_phi] + return k + + +def make_problem(): + nodes = [0.0, 1.0, 2.0] + meshgen = OrthogonalMeshGenerator(node_sets=[nodes, nodes, nodes]) + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + grid.SetOrthogonalBoundaries() + + xs = MultiGroupXS() + xs.LoadFromOpenSn("simple_5g_upscatter_fissile.xs") + + quad = GLCProductQuadrature3DXYZ(n_polar=2, n_azimuthal=4, scattering_order=0) + return DiscreteOrdinatesProblem( + mesh=grid, + num_groups=5, + groupsets=[ + { + "groups_from_to": [0, 0], + "angular_quadrature": quad, + "angle_aggregation_type": "single", + "angle_aggregation_num_subsets": 1, + "inner_linear_method": "petsc_gmres", + "l_abs_tol": 1.0e-10, + "l_max_its": 200, + "gmres_restart_interval": 30, + }, + { + "groups_from_to": [1, 1], + "angular_quadrature": quad, + "angle_aggregation_type": "single", + "angle_aggregation_num_subsets": 1, + "inner_linear_method": "petsc_gmres", + "l_abs_tol": 1.0e-10, + "l_max_its": 200, + "gmres_restart_interval": 30, + }, + { + "groups_from_to": [2, 2], + "angular_quadrature": quad, + "angle_aggregation_type": "single", + "angle_aggregation_num_subsets": 1, + "inner_linear_method": "petsc_gmres", + "l_abs_tol": 1.0e-10, + "l_max_its": 200, + "gmres_restart_interval": 30, + }, + { + "groups_from_to": [3, 3], + "angular_quadrature": quad, + "angle_aggregation_type": "single", + "angle_aggregation_num_subsets": 1, + "inner_linear_method": "petsc_gmres", + "l_abs_tol": 1.0e-10, + "l_max_its": 200, + "gmres_restart_interval": 30, + }, + { + "groups_from_to": [4, 4], + "angular_quadrature": quad, + "angle_aggregation_type": "single", + "angle_aggregation_num_subsets": 1, + "inner_linear_method": "petsc_gmres", + "l_abs_tol": 1.0e-10, + "l_max_its": 200, + "gmres_restart_interval": 30, + }, + ], + xs_map=[ + {"block_ids": [0], "xs": xs}, + ], + boundary_conditions=[ + {"name": "xmin", "type": "reflecting"}, + {"name": "xmax", "type": "reflecting"}, + {"name": "ymin", "type": "reflecting"}, + {"name": "ymax", "type": "reflecting"}, + {"name": "zmin", "type": "reflecting"}, + {"name": "zmax", "type": "reflecting"}, + ], + options={ + "use_precursors": False, + "verbose_inner_iterations": False, + "verbose_ags_iterations": False, + "verbose_outer_iterations": True, + "max_ags_iterations": 200, + "ags_tolerance": 1.0e-10, + }, + ) + + +def run_solver(name): + problem = make_problem() + if name == "PowerIteration": + solver = PowerIterationKEigenSolver(problem=problem, max_iters=200, k_tol=1.0e-9) + elif name == "SLEPcKrylovSchur": + solver = SLEPcKEigenSolver( + lbs_problem=problem, + max_iters=100, + eigen_tol=1.0e-9, + eps_type="krylovschur", + ) + elif name == "SLEPcPower": + solver = SLEPcKEigenSolver( + lbs_problem=problem, + max_iters=200, + eigen_tol=1.0e-9, + eps_type="power", + ) + else: + raise ValueError(f"Unknown solver {name}.") + + solver.Initialize() + solver.Execute() + k = solver.GetEigenvalue() + error = abs(k - EXPECTED_K) + if rank == 0: + print(f"{name} k-eigenvalue: {k:.12f}") + print(f"{name} k-error: {error:.6e}") + return name, k, error + + +if __name__ == "__main__": + num_procs = 4 + if size != num_procs: + sys.exit(f"Incorrect number of processors. Expected {num_procs} but got {size}.") + + reference_k = compute_reference_k() + if abs(reference_k - EXPECTED_K) > 1.0e-12: + raise RuntimeError( + f"Internal reference mismatch: computed {reference_k:.16e}, expected {EXPECTED_K:.16e}." + ) + + if rank == 0: + print(f"Expected k-eigenvalue: {EXPECTED_K:.12f}") + + failures = [] + for solver_name in ("PowerIteration", "SLEPcKrylovSchur", "SLEPcPower"): + name, k, error = run_solver(solver_name) + if error > 5.0e-6: + failures.append(f"{name} k={k:.12f} differs from expected {EXPECTED_K:.12f}.") + + if failures: + raise RuntimeError(" ; ".join(failures)) diff --git a/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_slepc_U235.py b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_slepc_U235.py new file mode 100644 index 0000000000..165c8f6dcc --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_slepc_U235.py @@ -0,0 +1,98 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +3D 172G KEigenvalue Solver test using power iteration and OpenMC MGXS library +Test: Final k-eigenvalue: 2.2800213 +""" + +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.xs import MultiGroupXS + from pyopensn.aquad import GLCProductQuadrature3DXYZ + from pyopensn.solver import DiscreteOrdinatesProblem, SLEPcKEigenSolver + +if __name__ == "__main__": + + num_procs = 4 + if size != num_procs: + sys.exit(f"Incorrect number of processors. Expected {num_procs} but got {size}.") + + # Mesh setup + Nx = 5 + Lx = 2.0 + xmin = 0.0 + dx = Lx / Nx + xmesh = [xmin + k * dx for k in range(Nx + 1)] + + Ny = 5 + Ly = 2.0 + ymin = 0.0 + dy = Ly / Ny + ymesh = [ymin + k * dy for k in range(Ny + 1)] + + Nz = 5 + Lz = 2.0 + zmin = 0.0 + dz = Lz / Nz + zmesh = [zmin + k * dz for k in range(Nz + 1)] + + meshgen = OrthogonalMeshGenerator(node_sets=[xmesh, ymesh, zmesh]) + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + + # Load cross section data using the OpenMC MGXS library + num_groups = 172 + xs_u235 = MultiGroupXS() + xs_u235.LoadFromOpenMC("u235_172g.h5", "u235", 294.0) + + # Solver Setup + phys = DiscreteOrdinatesProblem( + mesh=grid, + num_groups=num_groups, + groupsets=[ + { + "groups_from_to": (0, num_groups - 1), + "angular_quadrature": GLCProductQuadrature3DXYZ(n_polar=2, + n_azimuthal=4, + scattering_order=0), + "inner_linear_method": "petsc_gmres", + "l_max_its": 500, + "l_abs_tol": 1.0e-6, + }, + ], + xs_map=[ + {"block_ids": [0], "xs": xs_u235}, + ], + boundary_conditions=[ + {"name": "xmin", "type": "reflecting"}, + {"name": "xmax", "type": "reflecting"}, + {"name": "ymin", "type": "reflecting"}, + {"name": "ymax", "type": "reflecting"}, + {"name": "zmin", "type": "reflecting"}, + {"name": "zmax", "type": "reflecting"}, + ], + options={ + "use_precursors": False, + "verbose_inner_iterations": False, + "verbose_outer_iterations": True, + } + ) + k_solver = SLEPcKEigenSolver( + lbs_problem=phys, + max_iters=100, + eigen_tol=1.0e-8, + ) + k_solver.Initialize() + k_solver.Execute() + + k = k_solver.GetEigenvalue() + if rank == 0: + print(f"Python k-eigenvalue: {k}") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_slepc_power_U235.py b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_slepc_power_U235.py new file mode 100644 index 0000000000..9bb08d8d0b --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_slepc_power_U235.py @@ -0,0 +1,98 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +3D 172G KEigenvalue Solver test using power iteration and OpenMC MGXS library +Test: Final k-eigenvalue: 2.2800213 +""" + +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.xs import MultiGroupXS + from pyopensn.aquad import GLCProductQuadrature3DXYZ + from pyopensn.solver import DiscreteOrdinatesProblem, SLEPcKEigenSolver + +if __name__ == "__main__": + + num_procs = 4 + if size != num_procs: + sys.exit(f"Incorrect number of processors. Expected {num_procs} but got {size}.") + + # Mesh setup + Nx = 5 + Lx = 2.0 + xmin = 0.0 + dx = Lx / Nx + xmesh = [xmin + k * dx for k in range(Nx + 1)] + + Ny = 5 + Ly = 2.0 + ymin = 0.0 + dy = Ly / Ny + ymesh = [ymin + k * dy for k in range(Ny + 1)] + + Nz = 5 + Lz = 2.0 + zmin = 0.0 + dz = Lz / Nz + zmesh = [zmin + k * dz for k in range(Nz + 1)] + + meshgen = OrthogonalMeshGenerator(node_sets=[xmesh, ymesh, zmesh]) + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + + # Load cross section data using the OpenMC MGXS library + num_groups = 172 + xs_u235 = MultiGroupXS() + xs_u235.LoadFromOpenMC("u235_172g.h5", "u235", 294.0) + + # Solver Setup + phys = DiscreteOrdinatesProblem( + mesh=grid, + num_groups=num_groups, + groupsets=[ + { + "groups_from_to": (0, num_groups - 1), + "angular_quadrature": GLCProductQuadrature3DXYZ(n_polar=2, + n_azimuthal=4, + scattering_order=0), + "inner_linear_method": "petsc_gmres", + "l_max_its": 500, + "l_abs_tol": 1.0e-6, + }, + ], + xs_map=[ + {"block_ids": [0], "xs": xs_u235}, + ], + boundary_conditions=[ + {"name": "xmin", "type": "reflecting"}, + {"name": "xmax", "type": "reflecting"}, + {"name": "ymin", "type": "reflecting"}, + {"name": "ymax", "type": "reflecting"}, + {"name": "zmin", "type": "reflecting"}, + {"name": "zmax", "type": "reflecting"}, + ], + options={ + "use_precursors": False, + "verbose_inner_iterations": False, + "verbose_outer_iterations": True, + } + ) + k_solver = SLEPcKEigenSolver( + lbs_problem=phys, + max_iters=5000, + eigen_tol=1.0e-8, + eps_type="power" + ) + k_solver.Initialize() + k_solver.Execute() + + k = k_solver.GetEigenvalue() + print(f"Python k-eigenvalue: {k}") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_keigen/simple_5g_upscatter_fissile.xs b/test/python/modules/linear_boltzmann_solvers/transport_keigen/simple_5g_upscatter_fissile.xs index 13985e5213..9968a05848 100644 --- a/test/python/modules/linear_boltzmann_solvers/transport_keigen/simple_5g_upscatter_fissile.xs +++ b/test/python/modules/linear_boltzmann_solvers/transport_keigen/simple_5g_upscatter_fissile.xs @@ -1,5 +1,9 @@ +# Five-group fissile infinite-medium test data. +# The regression input uses five groupsets, one group per groupset. +# Any off-diagonal entry with gfrom > gto is groupset-to-groupset upscatter. NUM_GROUPS 5 NUM_MOMENTS 1 +NUM_PRECURSORS 0 SIGMA_T_BEGIN 0 1.0 @@ -46,29 +50,30 @@ GPRIME_G_VAL 4 4 0.002 PRODUCTION_MATRIX_END TRANSFER_MOMENTS_BEGIN +# Zeroth moment (l=0) M_GFROM_GTO_VAL 0 0 0 0.20 -M_GFROM_GTO_VAL 0 0 1 0.02 -M_GFROM_GTO_VAL 0 0 2 0.03 -M_GFROM_GTO_VAL 0 0 3 0.01 -M_GFROM_GTO_VAL 0 0 4 0.02 -M_GFROM_GTO_VAL 0 1 0 0.04 +M_GFROM_GTO_VAL 0 1 0 0.02 +M_GFROM_GTO_VAL 0 2 0 0.03 +M_GFROM_GTO_VAL 0 3 0 0.01 +M_GFROM_GTO_VAL 0 4 0 0.02 +M_GFROM_GTO_VAL 0 0 1 0.04 M_GFROM_GTO_VAL 0 1 1 0.25 -M_GFROM_GTO_VAL 0 1 2 0.02 -M_GFROM_GTO_VAL 0 1 3 0.03 -M_GFROM_GTO_VAL 0 1 4 0.01 -M_GFROM_GTO_VAL 0 2 0 0.06 -M_GFROM_GTO_VAL 0 2 1 0.05 +M_GFROM_GTO_VAL 0 2 1 0.02 +M_GFROM_GTO_VAL 0 3 1 0.03 +M_GFROM_GTO_VAL 0 4 1 0.01 +M_GFROM_GTO_VAL 0 0 2 0.06 +M_GFROM_GTO_VAL 0 1 2 0.05 M_GFROM_GTO_VAL 0 2 2 0.22 -M_GFROM_GTO_VAL 0 2 3 0.04 -M_GFROM_GTO_VAL 0 2 4 0.02 -M_GFROM_GTO_VAL 0 3 0 0.02 -M_GFROM_GTO_VAL 0 3 1 0.06 -M_GFROM_GTO_VAL 0 3 2 0.07 +M_GFROM_GTO_VAL 0 3 2 0.04 +M_GFROM_GTO_VAL 0 4 2 0.02 +M_GFROM_GTO_VAL 0 0 3 0.02 +M_GFROM_GTO_VAL 0 1 3 0.06 +M_GFROM_GTO_VAL 0 2 3 0.07 M_GFROM_GTO_VAL 0 3 3 0.24 -M_GFROM_GTO_VAL 0 3 4 0.03 -M_GFROM_GTO_VAL 0 4 0 0.01 -M_GFROM_GTO_VAL 0 4 1 0.03 -M_GFROM_GTO_VAL 0 4 2 0.04 -M_GFROM_GTO_VAL 0 4 3 0.05 +M_GFROM_GTO_VAL 0 4 3 0.03 +M_GFROM_GTO_VAL 0 0 4 0.01 +M_GFROM_GTO_VAL 0 1 4 0.03 +M_GFROM_GTO_VAL 0 2 4 0.04 +M_GFROM_GTO_VAL 0 3 4 0.05 M_GFROM_GTO_VAL 0 4 4 0.21 TRANSFER_MOMENTS_END diff --git a/test/python/modules/linear_boltzmann_solvers/transport_keigen/tests.json b/test/python/modules/linear_boltzmann_solvers/transport_keigen/tests.json index 0b8a8ec982..a4aa92c350 100644 --- a/test/python/modules/linear_boltzmann_solvers/transport_keigen/tests.json +++ b/test/python/modules/linear_boltzmann_solvers/transport_keigen/tests.json @@ -256,9 +256,64 @@ } ] }, + { + "file": "keigenvalue_transport_3d_slepc_U235.py", + "comment": + "3D 172G keigenvalue test using OpenMX MGXS cross-sections and SLEPc Krylov-Schur iteration", + "num_procs": 4, + "checks": [ + { + "type": "KeyValuePair", + "key": "Python k-eigenvalue:", + "goldvalue": 2.2800213, + "abs_tol": 1e-05 + } + ] + }, + { + "file": "keigenvalue_transport_3d_slepc_power_U235.py", + "comment": + "3D 172G keigenvalue test using OpenMX MGXS cross-sections and SLEPc power iteration", + "num_procs": 4, + "checks": [ + { + "type": "KeyValuePair", + "key": "Python k-eigenvalue:", + "goldvalue": 2.2800213, + "abs_tol": 1e-05 + } + ] + }, + { + "file": "keigenvalue_transport_3d_5g_upscatter_slepc.py", + "comment": + "3D 5G infinite-medium k-eigenvalue test with five groupsets, upscatter, analytic reference, and SLEPc", + "num_procs": 4, + "checks": [ + { + "type": "KeyValuePair", + "key": "PowerIteration k-eigenvalue:", + "goldvalue": 0.655371548443, + "abs_tol": 5e-06 + }, + { + "type": "KeyValuePair", + "key": "SLEPcKrylovSchur k-eigenvalue:", + "goldvalue": 0.655371548443, + "abs_tol": 5e-06 + }, + { + "type": "KeyValuePair", + "key": "SLEPcPower k-eigenvalue:", + "goldvalue": 0.655371548443, + "abs_tol": 5e-06 + } + ] + }, { "file": "keigenvalue_transport_3d_openmcxs_U235_smm.py", - "comment": "3D 84G keigenvalue test using OpenMC MGXS cross-sections and power iteration with SMM", + "comment": + "3D 84G keigenvalue test using OpenMC MGXS cross-sections and power iteration with SMM", "num_procs": 4, "checks": [ { @@ -374,7 +429,8 @@ }, { "file": "keigenvalue_transport_3d_openmcxs_UO2_crichardson.py", - "comment": "3D 172G KEigenvalue::Solver test using Power Iteration with Classic Richardson and OpenMC MGXS library", + "comment": + "3D 172G KEigenvalue::Solver test using Power Iteration with Classic Richardson and OpenMC MGXS library", "num_procs": 4, "checks": [ { @@ -403,7 +459,8 @@ }, { "file": "keigenvalue_transport_3d_2g_exo_wedge.py", - "comment": "3D_2G Infinite medium testing wedge elements of EXODUSII; should agree with the hex model", + "comment": + "3D_2G Infinite medium testing wedge elements of EXODUSII; should agree with the hex model", "num_procs": 1, "checks": [ { @@ -517,7 +574,8 @@ }, { "file": "keigenvalue_transport_3d_2g_exo_hex.py", - "comment": "3D_2G Infinite medium testing hex elements of EXODUSII; should agree with the wedge model", + "comment": + "3D_2G Infinite medium testing hex elements of EXODUSII; should agree with the wedge model", "num_procs": 1, "checks": [ { From 9c7ad6e8932071c1b68a436dfd5be7aaf88195a3 Mon Sep 17 00:00:00 2001 From: "W. Daryl Hawkins" Date: Wed, 29 Apr 2026 09:16:14 -0500 Subject: [PATCH 2/2] Updating SLEPc routines to use iterating logging methods. --- .../slepc_linear_keigen_solver.cc | 112 +++++++++++++++--- .../slepc_linear_keigen_solver.h | 5 + .../iterative_methods/iteration_logging.h | 20 ++++ .../xs}/simple_5g_upscatter_fissile.xs | 0 ...envalue_transport_3d_5g_upscatter_slepc.py | 5 +- ...value_transport_3d_5g_upscatter_solvers.py | 8 +- .../keigenvalue_transport_3d_slepc_U235.py | 2 +- ...igenvalue_transport_3d_slepc_power_U235.py | 2 +- .../transport_keigen/tests.json | 4 +- 9 files changed, 133 insertions(+), 25 deletions(-) rename test/{python/modules/linear_boltzmann_solvers/transport_keigen => assets/xs}/simple_5g_upscatter_fissile.xs (100%) diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/slepc_linear_keigen_solver.cc b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/slepc_linear_keigen_solver.cc index f08952c6c1..07aad64a99 100644 --- a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/slepc_linear_keigen_solver.cc +++ b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/slepc_linear_keigen_solver.cc @@ -5,6 +5,7 @@ #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/ags_linear_solver.h" #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/wgs_context.h" #include "modules/linear_boltzmann_solvers/lbs_problem/compute/lbs_compute.h" +#include "modules/linear_boltzmann_solvers/lbs_problem/iterative_methods/iteration_logging.h" #include "modules/linear_boltzmann_solvers/lbs_problem/vecops/lbs_vecops.h" #include "framework/math/petsc_utils/petsc_utils.h" #include "framework/math/math.h" @@ -16,8 +17,10 @@ #include #include #include +#include #include #include +#include namespace opensn { @@ -158,6 +161,18 @@ ApplyLinearKEigenOperator(SLEPcLinearKEigenContext& ctx, Vec x, Vec y) OpenSnLogicalErrorIf(not ags_solver, "SLEPcLinearKEigenSolver: AGS solver is not available."); ags_solver->Solve(); + ctx.last_ags_solve = ags_solver->GetLastSolveSummary(); + ctx.last_wgs_solves.clear(); + ctx.last_wgs_solves.reserve(do_problem->GetNumWGSSolvers()); + for (size_t gsid = 0; gsid < do_problem->GetNumWGSSolvers(); ++gsid) + { + auto wgs_context = + std::dynamic_pointer_cast(do_problem->GetWGSSolver(gsid)->GetContext()); + OpenSnLogicalErrorIf(not wgs_context, "SLEPcLinearKEigenSolver: Cast to WGSContext failed."); + if (HasIterationStatus(wgs_context->last_solve)) + ctx.last_wgs_solves.emplace_back(wgs_context->last_solve); + } + LBSVecOps::SetGroupScopedPETScVecFromPrimarySTLvector( *do_problem, 0, do_problem->GetNumGroups() - 1, y, phi_new_local); @@ -177,45 +192,75 @@ ShellMult(Mat M, Vec x, Vec y) return ApplyLinearKEigenOperator(*ctx, x, y); } -} // namespace +IterationStatus +SLEPcStatusToIterationStatus(const PetscInt nconv, + const PetscInt iterations, + const int maximum_iterations, + const IterationStatus inner_status) +{ + if (inner_status == IterationStatus::FAILED or inner_status == IterationStatus::LIMIT) + return inner_status; + if (nconv > 0) + return IterationStatus::CONVERGED; + return iterations >= maximum_iterations ? IterationStatus::LIMIT : IterationStatus::NONE; +} -static PetscErrorCode +PetscErrorCode SLEPcLinearKEigenMonitor(EPS unused_eps, PetscInt its, PetscInt unused_nconv, PetscScalar* unused_eigr, // NOLINT(readability-non-const-parameter) PetscScalar* unused_eigi, // NOLINT(readability-non-const-parameter) - PetscReal* errest, + PetscReal* errest, // NOLINT(readability-non-const-parameter) PetscInt nest, void* ctx) { (void)unused_eps; - (void)unused_nconv; (void)unused_eigr; (void)unused_eigi; const bool has_error_estimate = nest > 0 and errest != nullptr; + auto* kctx = static_cast(ctx); if (ctx) { - auto* kctx = static_cast(ctx); if (has_error_estimate) kctx->last_eps_residual = errest[0]; + kctx->last_eps_solve.num_iterations = static_cast(its); + kctx->last_eps_solve.metric_name = "residual"; + kctx->last_eps_solve.metric_value = kctx->last_eps_residual; + } + + if (kctx and kctx->do_problem->GetOptions().verbose_outer_iterations) + { + IterationSummary eps_summary = kctx->last_eps_solve; + if (not has_error_estimate) + eps_summary.metric_name.clear(); + + const auto ags_status = HasIterationStatus(kctx->last_ags_solve) ? kctx->last_ags_solve.status + : IterationStatus::NONE; + const auto inner_status = + MostSevereIterationStatus(ags_status, MostSevereIterationStatus(kctx->last_wgs_solves)); + const auto iteration_status = SLEPcStatusToIterationStatus( + unused_nconv, its, std::numeric_limits::max(), inner_status); + + log.Log() << program_timer.GetTimeString() << " " + << FormatEigenIteration("SLEPc EPS", + static_cast(its), + eps_summary, + kctx->last_ags_solve, + kctx->last_wgs_solves, + iteration_status); } - std::stringstream iter_info; - iter_info << program_timer.GetTimeString() << " SLEPc EPS iteration = " << its; - if (has_error_estimate) - iter_info << ", residual = " << std::scientific << std::setprecision(6) << errest[0] - << std::defaultfloat; - else - iter_info << ", status = iterating"; - log.Log() << iter_info.str(); return 0; } +} // namespace + void SLEPcLinearKEigenSolver::SetMonitor() { + EPSMonitorSet(eps_, SLEPcLinearKEigenMonitor, context_ptr_.get(), nullptr); } void @@ -275,12 +320,26 @@ SLEPcLinearKEigenSolver::Solve() EPSSetTolerances(eps_, tolerance_options.residual_absolute, tolerance_options.maximum_iterations); EPSSetType(eps_, eps_type_.c_str()); EPSSetInitialSpace(eps_, 1, &x_); - EPSMonitorSet(eps_, SLEPcLinearKEigenMonitor, context_ptr_.get(), nullptr); EPSSolve(eps_); // Check for convergence PetscInt nconv = 0; EPSGetConverged(eps_, &nconv); + PetscInt iterations = 0; + EPSGetIterationNumber(eps_, &iterations); + + const auto ags_status = + HasIterationStatus(ctx->last_ags_solve) ? ctx->last_ags_solve.status : IterationStatus::NONE; + const auto inner_status = + MostSevereIterationStatus(ags_status, MostSevereIterationStatus(ctx->last_wgs_solves)); + ctx->last_eps_solve.num_iterations = static_cast(iterations); + ctx->last_eps_solve.status = SLEPcStatusToIterationStatus( + nconv, iterations, tolerance_options.maximum_iterations, inner_status); + ctx->last_eps_solve.detail = + nconv > 0 ? "eigenpairs_converged=" + std::to_string(nconv) : "no_eigenpairs_converged"; + ctx->last_eps_solve.metric_name = "residual"; + ctx->last_eps_solve.metric_value = ctx->last_eps_residual; + if (nconv == 0) { OpenSnLogicalError(program_timer.GetTimeString() + " SLEPc EPS: No eigenpairs converged."); @@ -298,6 +357,9 @@ SLEPcLinearKEigenSolver::PreSolveCallback() auto ctx = std::dynamic_pointer_cast(context_ptr_); ctx->last_eps_residual = 1.0; ctx->last_operator_residual = 1.0; + ctx->last_eps_solve = {}; + ctx->last_ags_solve = {}; + ctx->last_wgs_solves.clear(); } void @@ -345,6 +407,28 @@ SLEPcLinearKEigenSolver::PostSolveCallback() OpenSnLogicalErrorIf(not std::isfinite(ctx->last_operator_residual), "SLEPcLinearKEigenSolver: final operator residual is not finite."); + if (do_problem->GetOptions().verbose_outer_iterations) + { + const auto ags_status = + HasIterationStatus(ctx->last_ags_solve) ? ctx->last_ags_solve.status : IterationStatus::NONE; + const auto inner_status = + MostSevereIterationStatus(ags_status, MostSevereIterationStatus(ctx->last_wgs_solves)); + if (inner_status == IterationStatus::FAILED or inner_status == IterationStatus::LIMIT) + ctx->last_eps_solve.status = inner_status; + + ctx->last_eps_solve.metric_name = "operator_residual"; + ctx->last_eps_solve.metric_value = ctx->last_operator_residual; + log.Log() << program_timer.GetTimeString() << " " + << FormatKEigenFinalSummary("SLEPc", + ctx->eigenvalue, + -1.0, + ctx->last_eps_solve.num_iterations, + "eps_iterations", + ctx->last_eps_solve.status); + log.Log() << program_timer.GetTimeString() << " " + << FormatIterationSummary("SLEPc EPS", ctx->last_eps_solve); + } + const auto ags_solver = do_problem->GetAGSSolver(); const double residual_tolerance = std::max({10.0 * tolerance_options.residual_absolute, ags_solver ? 10.0 * ags_solver->GetTolerance() : 0.0, diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/slepc_linear_keigen_solver.h b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/slepc_linear_keigen_solver.h index 91306002ad..5628a84d03 100644 --- a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/slepc_linear_keigen_solver.h +++ b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/slepc_linear_keigen_solver.h @@ -6,7 +6,9 @@ #include "framework/math/linear_solver/slepc_linear_eigen_solver.h" #include "framework/math/linear_solver/linear_solver_context.h" #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.h" +#include "modules/linear_boltzmann_solvers/lbs_problem/iterative_methods/iteration_logging.h" #include +#include namespace opensn { @@ -17,6 +19,9 @@ struct SLEPcLinearKEigenContext : public LinearEigenContext std::shared_ptr do_problem; double last_eps_residual = 1.0; double last_operator_residual = 1.0; + IterationSummary last_eps_solve; + IterationSummary last_ags_solve; + std::vector last_wgs_solves; explicit SLEPcLinearKEigenContext(std::shared_ptr do_problem) : do_problem(std::move(do_problem)) diff --git a/modules/linear_boltzmann_solvers/lbs_problem/iterative_methods/iteration_logging.h b/modules/linear_boltzmann_solvers/lbs_problem/iterative_methods/iteration_logging.h index f0815e3234..e5f4824084 100644 --- a/modules/linear_boltzmann_solvers/lbs_problem/iterative_methods/iteration_logging.h +++ b/modules/linear_boltzmann_solvers/lbs_problem/iterative_methods/iteration_logging.h @@ -182,6 +182,26 @@ FormatIterationSummary(const std::string& label, const IterationSummary& summary return out.str(); } +inline std::string +FormatEigenIteration(const std::string& phase_name, + const unsigned int iteration, + const IterationSummary& summary, + const IterationSummary& ags_summary, + const std::vector& wgs_summaries, + const IterationStatus outer_status) +{ + std::ostringstream out; + out << phase_name << " iteration = " << iteration; + if (not summary.metric_name.empty()) + AppendNumericField(out, summary.metric_name, summary.metric_value, Scientific(6)); + if (outer_status != IterationStatus::NONE) + out << ", status = " << IterationStatusName(outer_status); + if (HasIterationStatus(ags_summary)) + out << FormatNestedStatusSummary("AGS", ags_summary.status); + out << FormatNestedStatusCounts("WGS", wgs_summaries); + return out.str(); +} + inline std::string FormatKEigenOuterIteration(const std::string& phase_name, const unsigned int iteration, diff --git a/test/python/modules/linear_boltzmann_solvers/transport_keigen/simple_5g_upscatter_fissile.xs b/test/assets/xs/simple_5g_upscatter_fissile.xs similarity index 100% rename from test/python/modules/linear_boltzmann_solvers/transport_keigen/simple_5g_upscatter_fissile.xs rename to test/assets/xs/simple_5g_upscatter_fissile.xs diff --git a/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_5g_upscatter_slepc.py b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_5g_upscatter_slepc.py index cf07e6a249..e76758f577 100644 --- a/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_5g_upscatter_slepc.py +++ b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_5g_upscatter_slepc.py @@ -6,7 +6,7 @@ The homogeneous fully reflecting problem has a spatially constant solution. The reference k is the dominant eigenvalue of inv(T - S) F for the XS data in -simple_5g_upscatter_fissile.xs. +test/assets/xs/simple_5g_upscatter_fissile.xs. """ import os @@ -101,7 +101,7 @@ def make_problem(): grid.SetOrthogonalBoundaries() xs = MultiGroupXS() - xs.LoadFromOpenSn("simple_5g_upscatter_fissile.xs") + xs.LoadFromOpenSn("../../../../assets/xs/simple_5g_upscatter_fissile.xs") quad = GLCProductQuadrature3DXYZ(n_polar=2, n_azimuthal=4, scattering_order=0) return DiscreteOrdinatesProblem( @@ -173,7 +173,6 @@ def make_problem(): options={ "use_precursors": False, "verbose_inner_iterations": False, - "verbose_ags_iterations": False, "verbose_outer_iterations": True, "max_ags_iterations": 200, "ags_tolerance": 1.0e-10, diff --git a/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_5g_upscatter_solvers.py b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_5g_upscatter_solvers.py index 7d0d99f9b4..f8d446e202 100644 --- a/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_5g_upscatter_solvers.py +++ b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_5g_upscatter_solvers.py @@ -6,7 +6,7 @@ The homogeneous fully reflecting problem has a spatially constant solution. The reference k is the dominant eigenvalue of inv(T - S) F for the XS data in -simple_5g_upscatter_fissile.xs. +test/assets/xs/simple_5g_upscatter_fissile.xs. """ import os @@ -43,7 +43,7 @@ [0.030, 0.036, 0.012, 0.006, 0.005], [0.010, 0.012, 0.004, 0.002, 0.002], ] -EXPECTED_K = 0.647082071109 +EXPECTED_K = 0.6553715484430904 def solve_dense(matrix, rhs): @@ -73,7 +73,7 @@ def solve_dense(matrix, rhs): def apply_reference_operator(phi): loss = [ - [(SIGMA_T[g] if g == gp else 0.0) - SCATTER[gp][g] for gp in range(5)] + [(SIGMA_T[g] if g == gp else 0.0) - SCATTER[g][gp] for gp in range(5)] for g in range(5) ] fission_source = [ @@ -102,7 +102,7 @@ def make_problem(): grid.SetOrthogonalBoundaries() xs = MultiGroupXS() - xs.LoadFromOpenSn("simple_5g_upscatter_fissile.xs") + xs.LoadFromOpenSn("../../../../assets/xs/simple_5g_upscatter_fissile.xs") quad = GLCProductQuadrature3DXYZ(n_polar=2, n_azimuthal=4, scattering_order=0) groupsets = [] diff --git a/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_slepc_U235.py b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_slepc_U235.py index 165c8f6dcc..fdcb7a1652 100644 --- a/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_slepc_U235.py +++ b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_slepc_U235.py @@ -51,7 +51,7 @@ # Load cross section data using the OpenMC MGXS library num_groups = 172 xs_u235 = MultiGroupXS() - xs_u235.LoadFromOpenMC("u235_172g.h5", "u235", 294.0) + xs_u235.LoadFromOpenMC("../../../../assets/xs/u235_172g.h5", "u235", 294.0) # Solver Setup phys = DiscreteOrdinatesProblem( diff --git a/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_slepc_power_U235.py b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_slepc_power_U235.py index 9bb08d8d0b..c325705eab 100644 --- a/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_slepc_power_U235.py +++ b/test/python/modules/linear_boltzmann_solvers/transport_keigen/keigenvalue_transport_3d_slepc_power_U235.py @@ -51,7 +51,7 @@ # Load cross section data using the OpenMC MGXS library num_groups = 172 xs_u235 = MultiGroupXS() - xs_u235.LoadFromOpenMC("u235_172g.h5", "u235", 294.0) + xs_u235.LoadFromOpenMC("../../../../assets/xs/u235_172g.h5", "u235", 294.0) # Solver Setup phys = DiscreteOrdinatesProblem( diff --git a/test/python/modules/linear_boltzmann_solvers/transport_keigen/tests.json b/test/python/modules/linear_boltzmann_solvers/transport_keigen/tests.json index a4aa92c350..bc00d764cb 100644 --- a/test/python/modules/linear_boltzmann_solvers/transport_keigen/tests.json +++ b/test/python/modules/linear_boltzmann_solvers/transport_keigen/tests.json @@ -914,7 +914,7 @@ { "type": "KeyValuePair", "key": "PowerIteration k-eigenvalue:", - "goldvalue": 0.647082071109, + "goldvalue": 0.655371548443, "abs_tol": 5e-06 }, { @@ -926,7 +926,7 @@ { "type": "KeyValuePair", "key": "NonLinearKEigen k-eigenvalue:", - "goldvalue": 0.647082071109, + "goldvalue": 0.655371548443, "abs_tol": 5e-06 } ]