Migrate depletion module to C++ part 1 of N

[eepeterson]

Description

This PR is the first of a handful of steps towards moving the majority of the depletion module to C++. There are a number of reasons why we might want to do that some of which are listed below (there are tradeoffs as well obviously).

• Obviate the need to use the Python multiprocessing module which can cause installation and runtime headaches particularly on HPC systems.
• Faster CRAM solves without requiring SuperLU or Eigen libraries as additional dependencies.
• Easier tight coupling with multiphysics apps that need transmutation using OpenMC.

The general outline of the migration plan is as follows:

1. Port CRAM solves to C++ (This PR).
2. Port depletion chain infrastructure to C++
3. Port reaction rate calculation and matrix formation from tally result infrastructure to C++
4. Port integration schemes to C++
5. Refactor python side to simply orchestrate/manage the depletion calculation

Specifically this PR replaces CRAM16 and CRAM48 callables with analogous functions through openmc.lib and the C API. To make this work we need to replace scipy.sparse.linalg.spsolve calls with our own minimal sparse solver specific to CRAM. I debated trying to pull in SuperLU or Eigen for the sparse solve, but after digging into the details of both libraries and what is needed for CRAM specifically, I opted to roll our own solver which has many benefits. In particular, because the IPF form of CRAM has complex values on the diagonal whose imaginary parts are larger than 1, the pivots are unconditionally well-behaved and we don't need any of the partial (or complete) pivoting algorithms employed by SuperLU or Eigen.

The main additions of this PR are the CSCPattern and CSCMatrix classes for the SymbolicLUFactorization of matrices and data storage of burnup matrix elements. Then the BatemanSolver abstract base class and concrete implementation via the IPFCramSolver are added in C++ as mirrors of the current functionality on the Python side. The C++ CRAM solver is made accessible through openmc.lib and it is wired in to take the place of the scipy based solver for the CRAM16 and CRAM48 functions. The numeric linear algebra work specific to CRAM is kept in bateman_solvers.cpp since the numeric_factorize_cram function doesn't build the matrix $A dt - \theta_\ell I$ explicitly and isn't intended to be a generic LU factorization routine but rather does the numeric factorization on the fly specific to the CRAM matrix structure. The symbolic factorization of $A dt - \theta_\ell I$ is done once and reused by each numeric factorization for all the poles of the cram solve. Then triangular_solve_lu ingests the SymbolicLUFactorization and NumericLUFactorization and performs the linear solve for each pole which then get accumulated to compute the final nuclide density from the initial one.

Fixes # (issue)

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[yrrepy]

May be interesting to compare this with:
https://djsn.dev/post/lessons-porting-c++-cram-python/

[paulromano]

@eepeterson Thanks a ton for the PR! I haven't had a chance to look through in full yet, but right off the bat it would be very useful if you could share if you've done any comparisons on the accuracy of CRAM solutions for realistic systems as well as performance testing versus the current develop branch.

[shimwell]

@eepeterson Thanks a ton for the PR! I haven't had a chance to look through in full yet, but right off the bat it would be very useful if you could share if you've done any comparisons on the accuracy of CRAM solutions for realistic systems as well as performance testing versus the current develop branch.

This repo has a nice set of benchmarks, perhaps the SS316 steel is a reasonably complex one
https://github.com/jbae11/openmc_activator

[eepeterson]

Description

This PR is the first of a handful of steps towards moving the majority of the depletion module to C++. There are a number of reasons why we might want to do that some of which are listed below (there are tradeoffs as well obviously).

• Obviate the need to use the Python multiprocessing module which can cause installation and runtime headaches particularly on HPC systems.
• Faster CRAM solves without requiring SuperLU or Eigen libraries as additional dependencies.
• Easier tight coupling with multiphysics apps that need transmutation using OpenMC.

The general outline of the migration plan is as follows:

1. Port CRAM solves to C++ (This PR).
2. Port depletion chain infrastructure to C++
3. Port reaction rate calculation and matrix formation from tally result infrastructure to C++
4. Port integration schemes to C++
5. Refactor python side to simply orchestrate/manage the depletion calculation

Specifically this PR replaces CRAM16 and CRAM48 callables with analogous functions through openmc.lib and the C API. To make this work we need to replace scipy.sparse.linalg.spsolve calls with our own minimal sparse solver specific to CRAM. I debated trying to pull in SuperLU or Eigen for the sparse solve, but after digging into the details of both libraries and what is needed for CRAM specifically, I opted to roll our own solver which has many benefits. In particular, because the IPF form of CRAM has complex values on the diagonal whose imaginary parts are larger than 1, the pivots are unconditionally well-behaved and we don't need any of the partial (or complete) pivoting algorithms employed by SuperLU or Eigen.

The main additions of this PR are the CSCPattern and CSCMatrix classes for the SymbolicLUFactorization of matrices and data storage of burnup matrix elements. Then the BatemanSolver abstract base class and concrete implementation via the IPFCramSolver are added in C++ as mirrors of the current functionality on the Python side. The C++ CRAM solver is made accessible through openmc.lib and it is wired in to take the place of the scipy based solver for the CRAM16 and CRAM48 functions. The numeric linear algebra work specific to CRAM is kept in bateman_solvers.cpp since the numeric_factorize_cram function doesn't build the matrix $A dt - \theta_\ell I$ explicitly and isn't intended to be a generic LU factorization routine but rather does the numeric factorization on the fly specific to the CRAM matrix structure. The symbolic factorization of $A dt - \theta_\ell I$ is done once and reused by each numeric factorization for all the poles of the cram solve. Then triangular_solve_lu ingests the SymbolicLUFactorization and NumericLUFactorization and performs the linear solve for each pole which then get accumulated to compute the final nuclide density from the initial one.

Fixes # (issue)

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[yrrepy]

May be interesting to compare this with:
https://djsn.dev/post/lessons-porting-c++-cram-python/

[paulromano]

@eepeterson Thanks a ton for the PR! I haven't had a chance to look through in full yet, but right off the bat it would be very useful if you could share if you've done any comparisons on the accuracy of CRAM solutions for realistic systems as well as performance testing versus the current develop branch.

[shimwell]

@eepeterson Thanks a ton for the PR! I haven't had a chance to look through in full yet, but right off the bat it would be very useful if you could share if you've done any comparisons on the accuracy of CRAM solutions for realistic systems as well as performance testing versus the current develop branch.

This repo has a nice set of benchmarks, perhaps the SS316 steel is a reasonably complex one
https://github.com/jbae11/openmc_activator

[eepeterson]

Description

This PR is the first of a handful of steps towards moving the majority of the depletion module to C++. There are a number of reasons why we might want to do that some of which are listed below (there are tradeoffs as well obviously).

• Obviate the need to use the Python multiprocessing module which can cause installation and runtime headaches particularly on HPC systems.
• Faster CRAM solves without requiring SuperLU or Eigen libraries as additional dependencies.
• Easier tight coupling with multiphysics apps that need transmutation using OpenMC.

The general outline of the migration plan is as follows:

1. Port CRAM solves to C++ (This PR).
2. Port depletion chain infrastructure to C++
3. Port reaction rate calculation and matrix formation from tally result infrastructure to C++
4. Port integration schemes to C++
5. Refactor python side to simply orchestrate/manage the depletion calculation

Specifically this PR replaces CRAM16 and CRAM48 callables with analogous functions through openmc.lib and the C API. To make this work we need to replace scipy.sparse.linalg.spsolve calls with our own minimal sparse solver specific to CRAM. I debated trying to pull in SuperLU or Eigen for the sparse solve, but after digging into the details of both libraries and what is needed for CRAM specifically, I opted to roll our own solver which has many benefits. In particular, because the IPF form of CRAM has complex values on the diagonal whose imaginary parts are larger than 1, the pivots are unconditionally well-behaved and we don't need any of the partial (or complete) pivoting algorithms employed by SuperLU or Eigen.

The main additions of this PR are the CSCPattern and CSCMatrix classes for the SymbolicLUFactorization of matrices and data storage of burnup matrix elements. Then the BatemanSolver abstract base class and concrete implementation via the IPFCramSolver are added in C++ as mirrors of the current functionality on the Python side. The C++ CRAM solver is made accessible through openmc.lib and it is wired in to take the place of the scipy based solver for the CRAM16 and CRAM48 functions. The numeric linear algebra work specific to CRAM is kept in bateman_solvers.cpp since the numeric_factorize_cram function doesn't build the matrix $A dt - \theta_\ell I$ explicitly and isn't intended to be a generic LU factorization routine but rather does the numeric factorization on the fly specific to the CRAM matrix structure. The symbolic factorization of $A dt - \theta_\ell I$ is done once and reused by each numeric factorization for all the poles of the cram solve. Then triangular_solve_lu ingests the SymbolicLUFactorization and NumericLUFactorization and performs the linear solve for each pole which then get accumulated to compute the final nuclide density from the initial one.

Fixes # (issue)

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[yrrepy]

May be interesting to compare this with:
https://djsn.dev/post/lessons-porting-c++-cram-python/

[paulromano]

@eepeterson Thanks a ton for the PR! I haven't had a chance to look through in full yet, but right off the bat it would be very useful if you could share if you've done any comparisons on the accuracy of CRAM solutions for realistic systems as well as performance testing versus the current develop branch.

[shimwell]

@eepeterson Thanks a ton for the PR! I haven't had a chance to look through in full yet, but right off the bat it would be very useful if you could share if you've done any comparisons on the accuracy of CRAM solutions for realistic systems as well as performance testing versus the current develop branch.

This repo has a nice set of benchmarks, perhaps the SS316 steel is a reasonably complex one
https://github.com/jbae11/openmc_activator

[eepeterson]

Description

This PR is the first of a handful of steps towards moving the majority of the depletion module to C++. There are a number of reasons why we might want to do that some of which are listed below (there are tradeoffs as well obviously).

• Obviate the need to use the Python multiprocessing module which can cause installation and runtime headaches particularly on HPC systems.
• Faster CRAM solves without requiring SuperLU or Eigen libraries as additional dependencies.
• Easier tight coupling with multiphysics apps that need transmutation using OpenMC.

The general outline of the migration plan is as follows:

1. Port CRAM solves to C++ (This PR).
2. Port depletion chain infrastructure to C++
3. Port reaction rate calculation and matrix formation from tally result infrastructure to C++
4. Port integration schemes to C++
5. Refactor python side to simply orchestrate/manage the depletion calculation

Specifically this PR replaces CRAM16 and CRAM48 callables with analogous functions through openmc.lib and the C API. To make this work we need to replace scipy.sparse.linalg.spsolve calls with our own minimal sparse solver specific to CRAM. I debated trying to pull in SuperLU or Eigen for the sparse solve, but after digging into the details of both libraries and what is needed for CRAM specifically, I opted to roll our own solver which has many benefits. In particular, because the IPF form of CRAM has complex values on the diagonal whose imaginary parts are larger than 1, the pivots are unconditionally well-behaved and we don't need any of the partial (or complete) pivoting algorithms employed by SuperLU or Eigen.

The main additions of this PR are the CSCPattern and CSCMatrix classes for the SymbolicLUFactorization of matrices and data storage of burnup matrix elements. Then the BatemanSolver abstract base class and concrete implementation via the IPFCramSolver are added in C++ as mirrors of the current functionality on the Python side. The C++ CRAM solver is made accessible through openmc.lib and it is wired in to take the place of the scipy based solver for the CRAM16 and CRAM48 functions. The numeric linear algebra work specific to CRAM is kept in bateman_solvers.cpp since the numeric_factorize_cram function doesn't build the matrix $A dt - \theta_\ell I$ explicitly and isn't intended to be a generic LU factorization routine but rather does the numeric factorization on the fly specific to the CRAM matrix structure. The symbolic factorization of $A dt - \theta_\ell I$ is done once and reused by each numeric factorization for all the poles of the cram solve. Then triangular_solve_lu ingests the SymbolicLUFactorization and NumericLUFactorization and performs the linear solve for each pole which then get accumulated to compute the final nuclide density from the initial one.

Fixes # (issue)

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[yrrepy]

May be interesting to compare this with:
https://djsn.dev/post/lessons-porting-c++-cram-python/

[paulromano]

@eepeterson Thanks a ton for the PR! I haven't had a chance to look through in full yet, but right off the bat it would be very useful if you could share if you've done any comparisons on the accuracy of CRAM solutions for realistic systems as well as performance testing versus the current develop branch.

[shimwell]

@eepeterson Thanks a ton for the PR! I haven't had a chance to look through in full yet, but right off the bat it would be very useful if you could share if you've done any comparisons on the accuracy of CRAM solutions for realistic systems as well as performance testing versus the current develop branch.

This repo has a nice set of benchmarks, perhaps the SS316 steel is a reasonably complex one
https://github.com/jbae11/openmc_activator

[eepeterson]

Description

This PR is the first of a handful of steps towards moving the majority of the depletion module to C++. There are a number of reasons why we might want to do that some of which are listed below (there are tradeoffs as well obviously).

• Obviate the need to use the Python multiprocessing module which can cause installation and runtime headaches particularly on HPC systems.
• Faster CRAM solves without requiring SuperLU or Eigen libraries as additional dependencies.
• Easier tight coupling with multiphysics apps that need transmutation using OpenMC.

The general outline of the migration plan is as follows:

1. Port CRAM solves to C++ (This PR).
2. Port depletion chain infrastructure to C++
3. Port reaction rate calculation and matrix formation from tally result infrastructure to C++
4. Port integration schemes to C++
5. Refactor python side to simply orchestrate/manage the depletion calculation

Specifically this PR replaces CRAM16 and CRAM48 callables with analogous functions through openmc.lib and the C API. To make this work we need to replace scipy.sparse.linalg.spsolve calls with our own minimal sparse solver specific to CRAM. I debated trying to pull in SuperLU or Eigen for the sparse solve, but after digging into the details of both libraries and what is needed for CRAM specifically, I opted to roll our own solver which has many benefits. In particular, because the IPF form of CRAM has complex values on the diagonal whose imaginary parts are larger than 1, the pivots are unconditionally well-behaved and we don't need any of the partial (or complete) pivoting algorithms employed by SuperLU or Eigen.

The main additions of this PR are the CSCPattern and CSCMatrix classes for the SymbolicLUFactorization of matrices and data storage of burnup matrix elements. Then the BatemanSolver abstract base class and concrete implementation via the IPFCramSolver are added in C++ as mirrors of the current functionality on the Python side. The C++ CRAM solver is made accessible through openmc.lib and it is wired in to take the place of the scipy based solver for the CRAM16 and CRAM48 functions. The numeric linear algebra work specific to CRAM is kept in bateman_solvers.cpp since the numeric_factorize_cram function doesn't build the matrix $A dt - \theta_\ell I$ explicitly and isn't intended to be a generic LU factorization routine but rather does the numeric factorization on the fly specific to the CRAM matrix structure. The symbolic factorization of $A dt - \theta_\ell I$ is done once and reused by each numeric factorization for all the poles of the cram solve. Then triangular_solve_lu ingests the SymbolicLUFactorization and NumericLUFactorization and performs the linear solve for each pole which then get accumulated to compute the final nuclide density from the initial one.

Fixes # (issue)

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[yrrepy]

May be interesting to compare this with:
https://djsn.dev/post/lessons-porting-c++-cram-python/

[paulromano]

@eepeterson Thanks a ton for the PR! I haven't had a chance to look through in full yet, but right off the bat it would be very useful if you could share if you've done any comparisons on the accuracy of CRAM solutions for realistic systems as well as performance testing versus the current develop branch.

[shimwell]

@eepeterson Thanks a ton for the PR! I haven't had a chance to look through in full yet, but right off the bat it would be very useful if you could share if you've done any comparisons on the accuracy of CRAM solutions for realistic systems as well as performance testing versus the current develop branch.

This repo has a nice set of benchmarks, perhaps the SS316 steel is a reasonably complex one
https://github.com/jbae11/openmc_activator

[eepeterson]

Description

This PR is the first of a handful of steps towards moving the majority of the depletion module to C++. There are a number of reasons why we might want to do that some of which are listed below (there are tradeoffs as well obviously).

• Obviate the need to use the Python multiprocessing module which can cause installation and runtime headaches particularly on HPC systems.
• Faster CRAM solves without requiring SuperLU or Eigen libraries as additional dependencies.
• Easier tight coupling with multiphysics apps that need transmutation using OpenMC.

The general outline of the migration plan is as follows:

1. Port CRAM solves to C++ (This PR).
2. Port depletion chain infrastructure to C++
3. Port reaction rate calculation and matrix formation from tally result infrastructure to C++
4. Port integration schemes to C++
5. Refactor python side to simply orchestrate/manage the depletion calculation

Specifically this PR replaces CRAM16 and CRAM48 callables with analogous functions through openmc.lib and the C API. To make this work we need to replace scipy.sparse.linalg.spsolve calls with our own minimal sparse solver specific to CRAM. I debated trying to pull in SuperLU or Eigen for the sparse solve, but after digging into the details of both libraries and what is needed for CRAM specifically, I opted to roll our own solver which has many benefits. In particular, because the IPF form of CRAM has complex values on the diagonal whose imaginary parts are larger than 1, the pivots are unconditionally well-behaved and we don't need any of the partial (or complete) pivoting algorithms employed by SuperLU or Eigen.

The main additions of this PR are the CSCPattern and CSCMatrix classes for the SymbolicLUFactorization of matrices and data storage of burnup matrix elements. Then the BatemanSolver abstract base class and concrete implementation via the IPFCramSolver are added in C++ as mirrors of the current functionality on the Python side. The C++ CRAM solver is made accessible through openmc.lib and it is wired in to take the place of the scipy based solver for the CRAM16 and CRAM48 functions. The numeric linear algebra work specific to CRAM is kept in bateman_solvers.cpp since the numeric_factorize_cram function doesn't build the matrix $A dt - \theta_\ell I$ explicitly and isn't intended to be a generic LU factorization routine but rather does the numeric factorization on the fly specific to the CRAM matrix structure. The symbolic factorization of $A dt - \theta_\ell I$ is done once and reused by each numeric factorization for all the poles of the cram solve. Then triangular_solve_lu ingests the SymbolicLUFactorization and NumericLUFactorization and performs the linear solve for each pole which then get accumulated to compute the final nuclide density from the initial one.

Fixes # (issue)

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[yrrepy]

May be interesting to compare this with:
https://djsn.dev/post/lessons-porting-c++-cram-python/

[paulromano]

@eepeterson Thanks a ton for the PR! I haven't had a chance to look through in full yet, but right off the bat it would be very useful if you could share if you've done any comparisons on the accuracy of CRAM solutions for realistic systems as well as performance testing versus the current develop branch.

[shimwell]

@eepeterson Thanks a ton for the PR! I haven't had a chance to look through in full yet, but right off the bat it would be very useful if you could share if you've done any comparisons on the accuracy of CRAM solutions for realistic systems as well as performance testing versus the current develop branch.

This repo has a nice set of benchmarks, perhaps the SS316 steel is a reasonably complex one
https://github.com/jbae11/openmc_activator

[eepeterson]

Description

This PR is the first of a handful of steps towards moving the majority of the depletion module to C++. There are a number of reasons why we might want to do that some of which are listed below (there are tradeoffs as well obviously).

• Obviate the need to use the Python multiprocessing module which can cause installation and runtime headaches particularly on HPC systems.
• Faster CRAM solves without requiring SuperLU or Eigen libraries as additional dependencies.
• Easier tight coupling with multiphysics apps that need transmutation using OpenMC.

The general outline of the migration plan is as follows:

1. Port CRAM solves to C++ (This PR).
2. Port depletion chain infrastructure to C++
3. Port reaction rate calculation and matrix formation from tally result infrastructure to C++
4. Port integration schemes to C++
5. Refactor python side to simply orchestrate/manage the depletion calculation

Specifically this PR replaces CRAM16 and CRAM48 callables with analogous functions through openmc.lib and the C API. To make this work we need to replace scipy.sparse.linalg.spsolve calls with our own minimal sparse solver specific to CRAM. I debated trying to pull in SuperLU or Eigen for the sparse solve, but after digging into the details of both libraries and what is needed for CRAM specifically, I opted to roll our own solver which has many benefits. In particular, because the IPF form of CRAM has complex values on the diagonal whose imaginary parts are larger than 1, the pivots are unconditionally well-behaved and we don't need any of the partial (or complete) pivoting algorithms employed by SuperLU or Eigen.

The main additions of this PR are the CSCPattern and CSCMatrix classes for the SymbolicLUFactorization of matrices and data storage of burnup matrix elements. Then the BatemanSolver abstract base class and concrete implementation via the IPFCramSolver are added in C++ as mirrors of the current functionality on the Python side. The C++ CRAM solver is made accessible through openmc.lib and it is wired in to take the place of the scipy based solver for the CRAM16 and CRAM48 functions. The numeric linear algebra work specific to CRAM is kept in bateman_solvers.cpp since the numeric_factorize_cram function doesn't build the matrix $A dt - \theta_\ell I$ explicitly and isn't intended to be a generic LU factorization routine but rather does the numeric factorization on the fly specific to the CRAM matrix structure. The symbolic factorization of $A dt - \theta_\ell I$ is done once and reused by each numeric factorization for all the poles of the cram solve. Then triangular_solve_lu ingests the SymbolicLUFactorization and NumericLUFactorization and performs the linear solve for each pole which then get accumulated to compute the final nuclide density from the initial one.

Fixes # (issue)

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[yrrepy]

May be interesting to compare this with:
https://djsn.dev/post/lessons-porting-c++-cram-python/

[paulromano]

@eepeterson Thanks a ton for the PR! I haven't had a chance to look through in full yet, but right off the bat it would be very useful if you could share if you've done any comparisons on the accuracy of CRAM solutions for realistic systems as well as performance testing versus the current develop branch.

[shimwell]

@eepeterson Thanks a ton for the PR! I haven't had a chance to look through in full yet, but right off the bat it would be very useful if you could share if you've done any comparisons on the accuracy of CRAM solutions for realistic systems as well as performance testing versus the current develop branch.

This repo has a nice set of benchmarks, perhaps the SS316 steel is a reasonably complex one
https://github.com/jbae11/openmc_activator