Random Ray Kinetic Simulation Mode #3702
Conversation
- add `timestep_parameters` dictionary and `kinetic_simulation` switch to settings.py - add `bd_order` and `time_derivative_method` to `random_ray` dictionary in settings.py - add `density_timeseries` to material.py - add corresponding variables and machinery in the C++ core - add unit tests
- add some missing precursor tally machinery in reaction.cpp and tally.cpp, and output.cpp - add testing harness for kinetic simulations - remove some unuesd timers - implement time stepping loop in random_ray_simulation.cpp and supporting funtions.
- PROPOGATION -> PROPAGATION - Write random ray data to HDF5 - Write some extra HDF5 data - Typo fixes
…lations Kinetic simulations based on an eigenvalue solve need to multiply the prompt and delayed fission sources by the computed keff to bake-in the steady state initial condition. The keff scaling introduced in PR 3595 complicates this by multiplying tallied flux by keff. This is fine, however this new scaling factor is also applied to the final quantities which serve as the final estimate of quantities for each timestep in kinetic simulation. This would effetively undo the keff normalization putting the system out of steady state. In practice this looks like a shap step change in the values of tallied quantities. To address this, this commit introduces control flow to ensure that the keff scaling is applied to tallies at each timestep, but does not apply to the internal quantitities stored in memory until the final timestep, to ensure the keff normalization still bakes in that steady state.
- is_initial_condition - current_timetstep
yardasol
force-pushed
the
tdrr-clean-refactor
branch
from
3d4dc94 to
859fbec
Compare
yardasol
force-pushed
the
tdrr-clean-refactor
branch
from
1d72739 to
2f1b27b
Compare
yardasol
force-pushed
the
tdrr-clean-refactor
branch
from
85dc2a4 to
89abd44
Compare
|
Something is wrong with diagonal stabilization in the kinetic simulation mode. The simulation segfaults on the first k-eff simulation. Converting back to draft until I've identified and fixed the issue |
|
It seems a code block in |
| source_regions_.scalar_flux_new(se) = 0.0; | ||
| } | ||
|
|
||
| if (settings::kinetic_simulation && settings::create_delayed_neutrons) { |
There was a problem hiding this comment.
I believe the issue you were running into here has to do with the indexing of double precursors_new(int di) (this indexes precursors_new_ on the backend, which has length source_regions_.n_source_regions() * ndgroups_). n_source_elements() returns source_regions_.n_source_regions() * negroups_, which is notably quite a bit larger and so you are accidentally overrunning the precursor array.
There was a problem hiding this comment.
Good catch! Not sure how I missed that, I suppose that's what I get for coding late into the night. This can be fixed by swapping out n_source_elements() with n_delay_elements()
There was a problem hiding this comment.
It seems this flew under the radar previously as I was using more delay groups than energy groups for the verification. The indexing went way over when utilizing a 70-group structure, causing this odd behavior. Thanks again @nuclearkevin for your sharp eyes.
2 participants
Description
This PR adds a kinetic simulation mode to OpenMC's random ray solver to allow simulation of kinetic transients (i.e. time-dependent transients without temperature/reactivity feedback)
Specifically, this PR adds the following:
kinetic_simulationtoggle setting.CHI, to complementCHI_PROMPTandCHI_DELAYEDPRECURSORSrandom_ray_pin_cell.I realize this is quite a large feature and may take a while to review. I'd appreciate any discussion on necessary theory/user guide additions to the doc pages. Most of what can be done with the base random ray solver can be done with the kinetic simulation mode (with the exception of an incorrect treatment of void regions), so I figured not much past basic explanation of the feature and how to use it was needed.
Verification
I utilized two models to verify the kinetic simulation machinery
The pin cell model served as a quick sanity check for the reasonability of the my implementation’s solution. A fuel pin cell is easy to model and does not take a prohibitively long time to run. It is also relatively simple to model a pin cell using point-reactor kinetics, which serve as a robust check of my approach.
The more expensive C5G7 model was chosen due to
I utilize the same transient for both the pin cell model and the C5G7 model. The transient is based on Exercise 3-1 of the C5G7-TD benchmark. Moderator density is assumed to be at its nominal values at t=0 s . The moderator experiences a linear density decrease until it reaches a minimal value,$\omega=0.95$ for Exercise 3-1, at $t=1 s$ . The moderator then experiences a linear density increase until it returns to its nominal value at $t = 2s$ . The C5G7-TD reference results continue the simulation for another 8 seconds with the moderator at its nominal value. The moderator density changes imposed in the Exercise 3-1 transient, as well as the other three sub-cases for Exercise 3, are shown in the figure below:
Pincell Results
I compared the kinetic Random Ray implementation against PyRK, an open-source point-kinetics code. To generate the proper feedback coefficient, I ran several k-eigenvalue simulations along the range of densities in the Exercise 3-1 transient. Below are the PyRK and OpenMC results plotted against each other:
There is a large deviation between the time-dependent Random Ray solution and the point kinetics solution. While the two are very relatively close for the first leg of the density transient, once the density begins to increase the scalar flux and precursor solutions diverge. We know that the radial flux solution for a cylindrical reactor is the zeroth-order Bessel function$J_0$ , which decreases in magnitude as the radial point of interest moves away from the center of the reactor. The point kinetics flattens out this spatial dependence to a singular value, so it is expected that point kinetics overestimates both the scalar flux and the precursor concentration, as seen in the above figures. Similar statements can be made about the precursor solution.
I also wanted to compare the two methods (Time Isotropic (TI) and Source Derivative Propagation (SDP)). These two methods should be equivalent for this simple test case. TI is faster but will perform poorly near highly absorbing regions. SDP makes a higher order approximation for a greater execution cost.
The SDP method produces a slightly lower total neutron flux at the start of the transient, on the order of$10^-3$ below the TI scalar flux , but begins to approach the TI scalar flux until the inflection point at $t=1 s$ when it becomes slightly larger than the SDP solution. At $t=2 s$ , the SDP solution begins approaching the TI solution again before ending slightly below it for the remainder of the simulation. These errors are small enough that it is reasonable to conclude that the SDP and TI methods are equivalent for small models.
C5G7 results
I compared the results of the C5G7-TD transient run with OpenMC with reference data from Shen et al. (2019) on the same model and transient using the MOC code MPACT. The active length, inactive length, and rays per batch come from recommended parameters given by Tramm for the C5G7 model. The active and inactive batches are based on numbers used by Kraus et al. (2025). These simulation parameters are visible in the table below:
The figures below compare the MPACT reference solution to the OpenMC implementation of TI and SDP. All three methods are virtually indistinguishable during the density decrease and increase. The absolute error reaches a negative maximum of -0.02 and a positive maximum of 0.02 at the end of the density ramp. During the flat part of the density transient ( t=2 s and onward), the normalized fission rates of the SDP and TI methods are \sim0.05 below the MPACT normalized fission rate. Notice the oscillations in the absolute error. This is due to the stochastic nature of time-dependent Random Ray. The shape of the error plot and oscillations are very similar to Kraus’ findings.
Run Time Analysis
While the TI and SDP methods are equvialent, the TI method is around$\sim 5$ % and $\sim 19$ % faster than the SDP method for the pin cell and C5G7 Models.
Issues/Planned improvements
tally_delay_tasksobject inFlatSourceDomain::convert_source_regions_to_talliescopies a lot of existing code to constructtally_tasks. This shared code should ideally be parameterized in a helper function. An even better solution would be generalization oftally_tasksto include delay groupsopenmc_run_random_ray()function is getting quite long and has a lot of repeated code. I have a plan to parameterize this but this is for a future PR.Checklist