<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 3/10/17
Trying to get a critical infinite serpent simulation, e.g. $k_{\infty}$ = 1
U235 = .418%
U238 = .8625%
k = 1.07238
msr2g_enrU
2/10/17
Serpent run yielded k_eff of 1.03
msr2g_part_U_single_cell
2/10/17
Serpent run yielded k$_{eff}$ of 1.44
msr2g_part_U_full_core
2/13/17
Serpent run yielded k$_{eff}$ of .79216 with
- D ~ 55 cm
- H = 162.56 cm
- vacuum BCs
Dimensions based on buckling size given on pg. 39, MSRE Design and Operations, part iii, nuclear analysis
k$_{eff}$ = 1.07
D = 74.93 cm
H = 198.12 cm
vacuum BCs
graphite as special moderator material
k$_{eff}$ = 1.0545
D = 74.93 cm
H = 198.12 cm
vacuum BCs
graphite as regular doppler broadened material
)
1144 K fuel
k$_{eff}$ = 1.04149
Fuel temperature coefficient of reactivity = -3.08e-5 [$\frac{\delta k}{k} / ^{\circ}F$]
Step1:  1144 K fuel and graphite
k_eff = 1.02315
Step2:  Summary
Serpent Fuel C calculation
$\alpha_f = -3.08\cdot 10^{-5} [\frac{\delta k}{k} / ^{\circ}F]$
$\alpha_g = -4.34\cdot 10^{-5} [\frac{\delta k}{k} / ^{\circ}F]$
$\alpha_{tot} = -7.43\cdot 10^{-5} [\frac{\delta k}{k} / ^{\circ}F]$
MSRE report
$\alpha_f = -3.28\cdot 10^{-5} [\frac{\delta k}{k} / ^{\circ}F]$
$\alpha_g = -3.68\cdot 10^{-5} [\frac{\delta k}{k} / ^{\circ}F]$
$\alpha_{tot} = -6.96\cdot 10^{-5} [\frac{\delta k}{k} / ^{\circ}F]$
% error
$\alpha_f \rightarrow 5.96%$
$\alpha_g \rightarrow 18.15%$
$\alpha_{tot} \rightarrow 6.79%$
Step3:  msre_homogeneous
2/13/17
Serpent run yielded k$_{eff}$ of 1.03
General Notes
2/13/17
pg. 39, MSRE Design and Operations, part iii, nuclear analysis
Step4:  2/27/17
My two benchmarks in trying to understand serpent buckling are msre_homogeneous_critical_question_mark and msre_homog_rad_57_b1. These should both have the exact same material definitions and reactor height and differ only in the radius of the reactor. This is in fact the case, with the former having a radius of 73.86 cm and the latter having a radius of 57.1544 cm. What did this mean for respective $k_{inf}$? Intuitively we would hope for the change to be relatively small, since $\sigma(\vec{r},E)$ is unchanged between the two simulations. However, we could get changes because of the change in the reactor flux. So
Step5:  Above value perfectly matches one-group nubar
Step6:  Above matches perfectly with Serpent output B1_BUCKLING. So how does this compare with the original calculation where the material buckling was actually higher and yet the reactor was still sub-critical??
Step7:  So to move the reactor to a state of criticality (from a sub-critical state), the medium was made less absorptive, less fissile, and more diffusive. More particularly, $k_{\infty}$ was increased and the medium was made more diffusive. Ok, that seems like it could work.
See sage_scratch right before the beginning of 2/24/17 notes for some notes on buckling results. I believe that the quite large discrepancy between predicted critical buckling in Serpent and the actual buckling values is entirely due to large discrepancies in diffusion coefficient values. Let's consider the results from msre_homog_rad_57_b1_res.m. We have
Step8:  All the above results are very close to the b1 calculation results for the 57 cm reactor. So my question
<ASSISTANT_TASK:>
Python Code:
k_nom = 1.0545
k_f_1144 = 1.04149
fuel_reactivity = (k_f_1144 - k_nom) / k_nom / 400
print(fuel_reactivity)
Explanation: 3/10/17
Trying to get a critical infinite serpent simulation, e.g. $k_{\infty}$ = 1
U235 = .418%
U238 = .8625%
k = 1.07238
msr2g_enrU
2/10/17
Serpent run yielded k_eff of 1.03
msr2g_part_U_single_cell
2/10/17
Serpent run yielded k$_{eff}$ of 1.44
msr2g_part_U_full_core
2/13/17
Serpent run yielded k$_{eff}$ of .79216 with
- D ~ 55 cm
- H = 162.56 cm
- vacuum BCs
Dimensions based on buckling size given on pg. 39, MSRE Design and Operations, part iii, nuclear analysis
k$_{eff}$ = 1.07
D = 74.93 cm
H = 198.12 cm
vacuum BCs
graphite as special moderator material
k$_{eff}$ = 1.0545
D = 74.93 cm
H = 198.12 cm
vacuum BCs
graphite as regular doppler broadened material
)
1144 K fuel
k$_{eff}$ = 1.04149
Fuel temperature coefficient of reactivity = -3.08e-5 [$\frac{\delta k}{k} / ^{\circ}F$]
End of explanation
k_f_g_1144 = 1.02315
total_reactivity = (k_f_g_1144 - k_nom) / k_nom / 400
print(total_reactivity)
graph_reactivity = total_reactivity - fuel_reactivity
print(graph_reactivity)
from numpy import abs
print(abs(total_reactivity + 6.96e-5) / 6.96e-5 * 100)
print(abs(fuel_reactivity + 3.28e-5) / 3.28e-5 * 100)
print(abs(graph_reactivity + 3.68e-5) / 3.68e-5 * 100)
Explanation: 1144 K fuel and graphite
k_eff = 1.02315
End of explanation
from numpy import exp
alpha_fuel_faeh = 1.18e-4
alpha_fuel_kel = 1.8 * alpha_fuel_faeh
alpha_graph_faeh = 1.0e-5
alpha_graph_kel = 1.8 * alpha_graph_faeh
rho0_fuel = 2.146
rho0_graph = 1.86
T0 = 922
rho1144_fuel = rho0_fuel * exp(-alpha_fuel_kel * (1144 - T0))
rho1144_graph = rho0_graph * exp(-alpha_graph_kel * (1144 - T0))
print(rho1144_fuel)
print(rho1144_graph)
Explanation: Summary
Serpent Fuel C calculation
$\alpha_f = -3.08\cdot 10^{-5} [\frac{\delta k}{k} / ^{\circ}F]$
$\alpha_g = -4.34\cdot 10^{-5} [\frac{\delta k}{k} / ^{\circ}F]$
$\alpha_{tot} = -7.43\cdot 10^{-5} [\frac{\delta k}{k} / ^{\circ}F]$
MSRE report
$\alpha_f = -3.28\cdot 10^{-5} [\frac{\delta k}{k} / ^{\circ}F]$
$\alpha_g = -3.68\cdot 10^{-5} [\frac{\delta k}{k} / ^{\circ}F]$
$\alpha_{tot} = -6.96\cdot 10^{-5} [\frac{\delta k}{k} / ^{\circ}F]$
% error
$\alpha_f \rightarrow 5.96%$
$\alpha_g \rightarrow 18.15%$
$\alpha_{tot} \rightarrow 6.79%$
End of explanation
print(100 * .00272 / 1.02316)
Explanation: msre_homogeneous
2/13/17
Serpent run yielded k$_{eff}$ of 1.03
General Notes
2/13/17
pg. 39, MSRE Design and Operations, part iii, nuclear analysis: For one-region model of reactivity effects of temperature, geometric buckling used was that of a cylinder, 59 in. in diameter by 78 in. high. Three conditions:
salt and graphite at 1200 F
salt at 1600 F, graphite at 1200 F
salt and graphite at 1600 F
Temperature Reactivity coefficients both negative for fuel and graphite
Most definite reactor geometry and material descriptions: pages 18 & 19 of MSRE Design and Op., part iii
On those pages is a 2D R-Z 20-region model; within each region the material composition is considered to be uniform
Critical concentrations of Uranium were calculated for three different salt compositions
End of explanation
tot_gen_rate = 1.02331
tot_fiss_rate = 4.2e-1
tot_gen_rate / tot_fiss_rate
Explanation: 2/27/17
My two benchmarks in trying to understand serpent buckling are msre_homogeneous_critical_question_mark and msre_homog_rad_57_b1. These should both have the exact same material definitions and reactor height and differ only in the radius of the reactor. This is in fact the case, with the former having a radius of 73.86 cm and the latter having a radius of 57.1544 cm. What did this mean for respective $k_{inf}$? Intuitively we would hope for the change to be relatively small, since $\sigma(\vec{r},E)$ is unchanged between the two simulations. However, we could get changes because of the change in the reactor flux. So:
msre_homogeneous_critical_question_mark: abs_gc_kinf = 1.54774
msre_homog_rad_57_b1: abs_gc_kinf = 1.51718
I would contend that this is a small enough change to be satisfactory. Moreover, as calculated in the sage_scratch notebook, the critical material buckling changes by even less:
msre_homogeneous_critical_question_mark: 1.32e-3
msre_homog_rad_57_b1: 1.31e-3
Below generation and fission rates taken from msre_homogeneous_critical_question_mark_res.m:
End of explanation
b1_remxs = 2.37419e-3
b1_nsf = 3.66701e-3
b1_diff = 9.85681e-1
b1_mat_buckl = (b1_nsf - b1_remxs) / b1_diff
b1_k_inf = b1_nsf / b1_remxs
print(b1_mat_buckl)
print(b1_k_inf)
Explanation: Above value perfectly matches one-group nubar: 2.43645. It's obvious that in Serpent, the total loss rate is normalized to one, and consequently the generation rate is exactly equal to $k_{eff}$.
Relative loss mechanisms differ between the two different sizes as one would hope. Leakage more important in smaller reactor (oh really?):
57 cm:
TOT_ABSRATE               (idx, [1:   2]) = [  5.40995E-01 0.00021 ];
TOT_LEAKRATE              (idx, [1:   2]) = [  4.59005E-01 0.00025 ];
74 cm:
TOT_ABSRATE               (idx, [1:   2]) = [  6.61174E-01 0.00017 ];
TOT_LEAKRATE              (idx, [1:   2]) = [  3.38826E-01 0.00033 ];
I don't even think that the idea of a critical buckling "search" makes any bloody sense. That implies to me that you are going to change the materials! To me the search should consist of setting $k_{eff}$ equal to one, and then calculating the appropriate geometric buckling. That's what makes sense to me. I just don't think that the buckling search works in Serpent. I decreased the geometric buckling, attempting to reach a value of 1.31e-3, starting from 3.52e-3 for the 57 cm radius reactor but didn't get there...I only got down to 1.90e-3 with the reactor radius of 74 cm. However, that was plenty far enough to go below the material buckling resulting in a super-critical reactor. That's why I don't think the buckling method in Serpent works.
Ok, we're at assembly level (or some other arbitrary level). Typically, even at assembly level we use reflective boundary conditions at first to isolate the assembly from the rest of the reactor core. Treating the current level this way yields a $k_{inf}$ of sorts and a spectrum that we can call an infinite medium spectrum because of the reflective conditions. However, we want to determine a criticality spectrum. A critical system can be realized in a few ways. One way would be to add some amount of absorber (if $k_{inf}>1$; negative absorber if the system is sub-critical). However, a more physically realistic way is to add leakage; e.g. if the reactor core is critical but one piece's $k_{inf}$ is subcritical, then there will be neat leakage into that component to balance.
Two group discretization results in two eigenvalues and their associated eigenvectors (e.g. two eigenvectors total) as long as the matrix has full rank.
It might be critical to understand that in general reactor theory language when someone talks about the neutron spectrum, s/he is talking about the neutron energy distribution and may very well not be referring at all to the spatial distribution of neutrons. Spectrum -> energy. And so an infinite spectrum may refer only to the neutron energy distribution produced by an infinite medium, and consequently the critical spectrum may only refer to the change in the neutron energy distribution from the infinite spectrum. E.g. we construct a relation between the leakage and infinite spectra (Stammler pg. 356):
\begin{equation}
r_L = r_{\infty} \frac{1+x}{1 + xk_{\infty}}
\end{equation}
where $r = \frac{\phi_1}{\phi_2}$, $r_{\infty} = \frac{\Sigma_{12}}{\Sigma_{a2}}$ and $x = L^2/\tau$ where:
\begin{equation}
L^2 = \frac{D_2}{\Sigma_{a2}}\
\tau = \frac{D_1}{\Sigma_{12} + \Sigma_{a1}}
\end{equation}
Note that the above equation for $r_L$ is valid only when the overall reactor is critical and only thermal neutrons cause fission events. A more general equation for the spectrum index is:
\begin{equation}
r_i = \frac{r_{\infty}}{1+L^2 B_i^2}
\end{equation}
where i equals either 1 or 2 for the two group case and the equations for $B_1^2$ and $B_2^2$ can be found on page 355 of Stammler. Spectrum eigenvalue problem corresponds to material buckling (largest of which is $B_1^2$ and is called the fundamental or asymptotic mode; for two group problems $B_2^2$ is referred to as the transient mode; presumably for a G group problem, anything other than $B_1^2$ would be referred to as a transient mode.) Spatial eigenvalue problem corresponds to geometric buckling. Here's $B_1^2$:
\begin{equation}
B_1^2 = \frac{\lambda k_{\infty} - 1}{\tau + \frac{\tau}{\tau*}L^2}
\end{equation}
With only thermal fission, $\tau* = \tau$. Note that $B_1^2$ is a function only of the material data and the multiplication factor of the reactor. In general, the corresponding Helmholtz equation for the spatial shape can be solved with appropriate boundary conditions determined by continuity of flux and current (assuming we are currently only examining one part of the reactor). However, if we require that the flux in all groups must vanish at the boundaries of a bare homogeneous system, then the Helmholtz equation for the spatial shape becomes an eigenvalue problem.
It should be noted from a math standpoint that eigenvalue/eigenfunction problems can only arise in differential equations when boundary conditions are homogeneous.
What I want to figure out tomorrow is exactly how the calculation of few group cross sections is implemented in serpent; e.g. look at the actual code calculation routines.
3/1/17
Homogenized cross section of reaction $i$:
\begin{equation}
i = j + 1
\end{equation}
\begin{equation}
\Sigma_{i,g} = \frac{\int_V \int_{E_g}^{E_{g-1}} \Sigma_i(\vec{r}, E)\ \phi(\vec{r},E)\ d^3r\ dE}{\int_V \int_{E_g}^{E_{g-1}} \phi(\vec{r},E)\ d^3r\ dE} \quad eq. 17.3.1.1
\end{equation}
Deterministic transport version:
\begin{equation}
\Sigma_{i,g} = \frac{ \sum_j \left( V_j \sum_h \left( \Sigma_{i,j,h} \Phi_{j,h} \right) \right) } 
{ \sum_j \left( V_j \sum_h \Phi_{j,h} \right) } \quad eq. 17.3.1.2
\end{equation}
In a cell with reflective boundary conditions and homogeneous material composition, the flux has no spatial dependence and is thus constant with respect to the spatial coordinates. However, it will not be constant with respect to a continuous energy variable. In general the determinstic transport calculation proceeds like this, similar to that outlined in Figure I.3 in Stammler (pg. 42):
Within homogeneous material regions, micro-group cross sections are calculated via Monte Carlo methods (which use coninuous energy cross sections). Double check that these micro-group cross sections are indeed calculated with Monte Carlo methods. Within each homogeneous region: $\Sigma_i(\vec{r},E) \rightarrow \Sigma_i(E)$.
Armed with micro-group cross sections in each region, the micro-group fluxes in a unit/pin-cell are calculated from some analytical or numerical transport method assuming reflective boundary conditions at the outer unit/pin-cell boundary (e.g. imagine the calculations perfomed in the 501 neutron project but with less approximations, e.g. we're definitely not doing diffusion). If a unit/pin-cell only has one material, then the micro-group flux solutions will display no spatial dependence.
Or put another way:
This step needs to be studied more. Monte Carlo calculation yields continuous energy neutron flux (restricted to what region?). This continuous energy flux is used to collapse continuous energy cross sections in homogeneous material regions to micro-group cross sections using continuos energy neutron flux
Unit/pin-cell calculation yields micro-group neutron fluxes within the unit/pin-cell, using reflective boundary conditions to completely isolate the unit/pin-cell from the outside world. If the unit/pin-cell is heterogeneous (e.g. multiple material regions), then the micro-group fluxes will display spatial dependence. If so, then a discretized form of equation 17.3.1.1 (using only a spatial filter) can be used to convert homogeneous micro-group cross sections into cell-averaged micro-group cross sections. Similarly, spatially-varying pin-cell fluxes can be homogenized by simple cell averaging (e.g. integrate cell flux over cell volume and then divide by cell volume) to create space indepenent cell average fluxes. Micro-group cell-averaged fluxes are then used to collapse cell-averaged micro-group cross sections into few-group, assembly-averaged cross sections using equation 17.3.1.2. (May also do an intermediate calculation with B2 groups).
Another very useful diagram is Stammler X.4. It's perhaps the best one I've seen
3/7/17
Understanding Serpent output:
DIFFAREA = RES_FG_L2
MAT_BUCKLING = RES_MAT_BUCKL
ABS_GC_KINF = RES_ABS_GC_KINF
Calculation of RES_MAT_BUCKL (where is the index looping over the few groups):
c
      if ((div = BufVal(RES_FG_L2, i)) &gt; 0.0)
        {
          val = (BufVal(RES_ABS_GC_KINF, i) - 1.0)/div;
          AddStat(RES_MAT_BUCKL, i, val);
Note that just because BufVal takes an index doesn't mean that the first argument is a true array, e.g. RES_ABS_GC_KINF only has one value (e.g. one mean and one measure of variance). 
Ok, here is how RES_ABS_GC_KINF is calculated:
c
      if ((loss = abs - n2n) &gt; 0.0)
    {
      keff = src/loss;
      AddStat(RES_ABS_GC_KINF, i, keff);
    }
End of explanation
remxs = 2.21079e-3
nsf = 3.35432e-3
diff = 5.31584e-1
mat_buckl = (nsf - remxs) / diff
k_inf = nsf / remxs
print(mat_buckl)
print(k_inf)
Explanation: Above matches perfectly with Serpent output B1_BUCKLING. So how does this compare with the original calculation where the material buckling was actually higher and yet the reactor was still sub-critical??
End of explanation
b1_73_remxs = 2.38026e-3
b1_73_nsf = 3.67837e-3
b1_73_diff = 9.85825e-1
mat_buckl_73 = (b1_73_nsf - b1_73_remxs) / b1_73_diff
k_inf_73 = b1_73_nsf / b1_73_remxs
print(mat_buckl_73)
print(k_inf_73)
Explanation: So to move the reactor to a state of criticality (from a sub-critical state), the medium was made less absorptive, less fissile, and more diffusive. More particularly, $k_{\infty}$ was increased and the medium was made more diffusive. Ok, that seems like it could work.
See sage_scratch right before the beginning of 2/24/17 notes for some notes on buckling results. I believe that the quite large discrepancy between predicted critical buckling in Serpent and the actual buckling values is entirely due to large discrepancies in diffusion coefficient values. Let's consider the results from msre_homog_rad_57_b1_res.m. We have:
DIFFCOEFF = 5.31584E-1
LEAK_DIFFCOEFF = 8.72571E-1
P1_DIFFCOEFF = 9.37823E-1
B1_DIFFCOEFF = 9.85681E-1
If we create a cylindrical reactor that from diffusion theory should have a geometric buckling equal to the predicted critical buckling predicted in the 57 cm radius Serpent simulation, then we produce this set of diffusion coefficients with Serpent:
DIFFCOEFF = 6.46163E-1
LEAK_DIFFCOEFF = 6.04551E-1
P1_DIFFCOEFF = 9.30344E-1
B1_DIFFCOEFF = 9.85825E-1
Here's also the important b1 results for the 73 cm radius (supposed to be critical but isn't) simulation:
End of explanation
40 * 20000
(381.99 + 160.60) / 4
381.99 + 160.60
Explanation: All the above results are very close to the b1 calculation results for the 57 cm reactor. So my question: how is the b1 diffusion coefficient calculated and how does it differ from calculation of DIFFCOEFF?
Serpent 1 and 2 results for the same input file are pretty much identical, which I suppose should be expected. However, the Serpent 2 result file does contain a lot more information.
3/9/17
Analog estimator: score the physical interactions (fission, capture, etc.) or event sequences that are simulated during the calculation. Analog estimates are directly related to the simulation process.
Note: in general when the chain of equalities ends, that generally means I've reached a point where all arguments are generated via Monte Carlo analog or implicit estimates, e.g. none of the arguments are derived quantities, e.g. if I search in collectresults.c, their first hit will be on the RHS of an equality sign. 
Default calculation of $D$:
\begin{equation}
D_g = \frac{\bar{r_g^2}}{6}\ \Sigma_{r,g}
\end{equation}
where $\bar{r_g^2}$ is the mean value of all scored square distances that neutrons have made within energy group $g$, which is an analog estimator. The starting points are the locations where neutrons have entered the group by fission or scattering from another (higher) energy group. The end points are the locations at which the neutrons are lost from the group either by absorption or scattering to another (lower) energy group.
Leakage calculation of $D$:
\begin{equation}
D_g = \frac{leak}{B_m^2}\
B_m^2 = \frac{k_{\infty} - 1}{L_g^2}\
L_g^2 = \frac{\bar{r_g^2}}{6}\
k_{\infty} = \frac{src}{loss} = \frac{\Sigma_{g,f}\bar{\nu}}{\Sigma_{g,f} + \Sigma_{g,c} - \Sigma_{n\rightarrow n}}
\end{equation}
where $leak$ comes directly from Monte Carlo scoring of neutrons exiting the simulation domain I believe.
P1 calculation of $D$:
\begin{equation}
D_g = \frac{1}{3\Sigma_{g,tr}}\
\Sigma_{g,tr} = \Sigma_{g,t} - \Sigma_{g,s1}
\end{equation}
B1 calculation of $D$ (note that this occurs in b1solver.c which is different from all the other calculations of $D$ that occur in collectresults.c; moreover the below equation comes from the Fridman and Lepannen PHYSOR paper as opposed to directly from the source code):
\begin{equation}
D_g = \frac{J_g}{|B|\phi_g}
\end{equation}
where $J_g$ is the group current and $|B|$ is the magnitude of the calculated critical buckling.
$\Sigma_r$
$\bar{\nu}\Sigma_f$
$D$
Ok, B1 calculations just don't make sense in multiple regions in Serpent. Serpent performs criticality searches in both regions independently. So that's out!
Scratch
End of explanation
<END_TASK>