

In [1]:

    
cd ~/projects/Okapi/tests/transfers/okapi-moose/pin_coupling/









    



/Users/lindad/projects/Okapi/tests/transfers/okapi-moose/pin_coupling



In [2]:

    
import openmc

import matplotlib.pyplot as plt
%matplotlib inline

Run with 40 active batches, 100k particles per batch



In [47]:

    
with openmc.StatePoint('statepoint.50.h5') as sp:
    mesh_tally = sp.get_tally(id=2)
    zernike_tally = sp.get_tally(id=1) 

scores = ['kappa-fission', 'fission', 'flux']
mesh_tallies = [mesh_tally.get_slice(scores=[score]) for score in scores]
for tally in mesh_tallies:
    tally.mean.shape = (100, 100)

for index, tally in enumerate(mesh_tallies):
    fig = plt.subplot(131 + index)
    fig.imshow(tally.mean)



In [48]:

    
import numpy as np

mesh_r = np.linspace(-60, 60, 100)

plt.plot(mesh_r, mesh_tallies[0].mean[50,:], label="kappa-fission")
plt.legend()









    Out[48]:





<matplotlib.legend.Legend at 0x1277db630>



In [49]:

    
plt.plot(mesh_r, mesh_tallies[1].mean[50,:], label="fission")
plt.legend()









    Out[49]:





<matplotlib.legend.Legend at 0x12751f550>



In [50]:

    
plt.plot(mesh_r, mesh_tallies[2].mean[50,:], label="flux")
plt.legend()









    Out[50]:





<matplotlib.legend.Legend at 0x1275d6c50>

Run with 10 active batches, 100k particles per batch



In [43]:

    
with openmc.StatePoint('statepoint.20.h5') as sp:
    mesh_tally = sp.get_tally(id=2)
    zernike_tally = sp.get_tally(id=1) 

scores = ['kappa-fission', 'fission', 'flux']
mesh_tallies = [mesh_tally.get_slice(scores=[score]) for score in scores]
for tally in mesh_tallies:
    tally.mean.shape = (100, 100)

for index, tally in enumerate(mesh_tallies):
    fig = plt.subplot(131 + index)
    fig.imshow(tally.mean)



In [44]:

    
import numpy as np

mesh_r = np.linspace(-60, 60, 100)

plt.plot(mesh_r, mesh_tallies[0].mean[50,:], label="kappa-fission")
plt.legend()









    Out[44]:





<matplotlib.legend.Legend at 0x126b09c18>



In [45]:

    
plt.plot(mesh_r, mesh_tallies[1].mean[50,:], label="fission")
plt.legend()









    Out[45]:





<matplotlib.legend.Legend at 0x126cfeeb8>



In [46]:

    
plt.plot(mesh_r, mesh_tallies[2].mean[50,:], label="flux")
plt.legend()









    Out[46]:





<matplotlib.legend.Legend at 0x1266cd400>

Run with 10 active batches, 10k particles per batch



In [59]:

    
with openmc.StatePoint('statepoint.20.h5') as sp:
    mesh_tally = sp.get_tally(id=2)
    zernike_tally = sp.get_tally(id=1) 

scores = ['kappa-fission', 'fission', 'flux', 'total']
mesh_tallies = [mesh_tally.get_slice(scores=[score]) for score in scores]
for tally in mesh_tallies:
    tally.mean.shape = (100, 100)

for index, tally in enumerate(mesh_tallies):
    fig = plt.subplot(101 + 10 * len(scores) + index)
    fig.imshow(tally.mean)



In [60]:

    
import numpy as np

mesh_r = np.linspace(-60, 60, 100)

plt.plot(mesh_r, mesh_tallies[0].mean[50,:], label="kappa-fission")
plt.legend()









    Out[60]:





<matplotlib.legend.Legend at 0x12740df60>



In [61]:

    
plt.plot(mesh_r, mesh_tallies[1].mean[50,:], label="fission")
plt.legend()









    Out[61]:





<matplotlib.legend.Legend at 0x127ecfd68>



In [62]:

    
plt.plot(mesh_r, mesh_tallies[2].mean[50,:], label="flux")
plt.legend()









    Out[62]:





<matplotlib.legend.Legend at 0x127f3fc18>



In [63]:

    
plt.plot(mesh_r, mesh_tallies[3].mean[50,:], label="total")
plt.legend()









    Out[63]:





<matplotlib.legend.Legend at 0x127f92240>

The above calculations show that the tallies really are normalized to the number of source neutrons. This is done through the total_weight variable in tally_header.F90



In [56]:

    
mesh_tallies[0].num_realizations









    Out[56]:





10



In [64]:

    
.01 * 10000









    Out[64]:





100.0

Zernike calculations



In [65]:

    
zernike_fission = zernike_tally.get_slice(scores=['kappa-fission'])
zernike_flux = zernike_tally.get_slice(scores=['flux'])



In [66]:

    
from series.Zernike import Zernike

max_order = 10
zernike_list = [Zernike(order, 0, 0, 60) for order in range(max_order + 1)]

zernike_flux.mean.shape = (zernike_flux.mean.size)

for order in range(max_order + 1):
    n = (order + 1) * (order + 2) // 2
    zernike_list[order].coefficients = zernike_flux.mean[:n]

import numpy as np

r = np.linspace(0, 60, 100)
z_evals = [np.zeros((r.size)) for zernike in zernike_list]

for r_index, radius in enumerate(r):
    for z_index, zernike in enumerate(zernike_list):
        z_evals[z_index][r_index] = zernike(radius, 0, normalization="orthonormal")

for order, z in enumerate(z_evals):
    if order % 2 == 0 and order <= 4:
        plt.plot(r, z, label=str(order))
plt.legend()









    Out[66]:





<matplotlib.legend.Legend at 0x126c98358>



In [68]:

    
zernike_list[0].coefficients









    Out[68]:





array([101.14264761])



In [71]:

    
from math import pi

z_evals[0][0] * pi









    Out[71]:





101.14264760764684



In [82]:

    
for z in z_evals:
    print(np.trapz(z * r, r) / np.trapz(r, r) * pi)









    



101.14264760764684
100.05844990007022
100.15523528161962
100.43465805412225
76.84860401899553
77.37545672633946
109.28574325268279
109.7084986885324
51.18621682319467
51.69695846172283
179.95786331423693

Above we see the role of uncertainty beginning to set it

The first Zernike coefficient should match a tally with no filters.



In [78]:

    
np.trapz(z_evals[0] * r, r) / np.trapz(r, r) * pi









    Out[78]:





101.14264760764684



In [79]:

    
np.trapz(z_evals[1] * r, r) / np.trapz(r, r) * pi









    Out[79]:





100.05844990007022



In [80]:

    
np.trapz(z_evals[2] * r, r) / np.trapz(r, r) * pi









    Out[80]:





100.15523528161962



In [81]:

    
np.trapz(z_evals[3] * r, r) / np.trapz(r, r) * pi









    Out[81]:





100.43465805412225

Below, note the oscillations introduced by a high order Zernike



In [54]:

    
from math import cos, pi

def sinusoid(r, phi, R=60):
    return cos(pi * r / (2 * R))

for zernike in zernike_list:
    zernike.generateCoefficients(sinusoid, num_rintervals=100, num_aintervals=10, normalization="standard")

for r_index, radius in enumerate(r):
    for z_index, zernike in enumerate(zernike_list):
        z_evals[z_index][r_index] = zernike(radius, 0, normalization="orthonormal")

for order, z in enumerate(z_evals):
    if order == 10 or order == 9:
        plt.plot(r, z, label=str(order))
true = np.zeros((r.size))
for index, radius in enumerate(r):
    true[index] = sinusoid(radius, 0)
plt.plot(r, true, label='true')
plt.legend()









    Out[54]:





<matplotlib.legend.Legend at 0x11e3504a8>

These oscillations can be removed by increasing the accuracy of the integration performed for coefficient generation



In [55]:

    
from math import cos, pi

def sinusoid(r, phi, R=60):
    return cos(pi * r / (2 * R))

for zernike in zernike_list:
    zernike.generateCoefficients(sinusoid, num_rintervals=100, num_aintervals=20, normalization="standard")

for r_index, radius in enumerate(r):
    for z_index, zernike in enumerate(zernike_list):
        z_evals[z_index][r_index] = zernike(radius, 0, normalization="orthonormal")

for order, z in enumerate(z_evals):
    if order == 10 or order == 9:
        plt.plot(r, z, label=str(order))
true = np.zeros((r.size))
for index, radius in enumerate(r):
    true[index] = sinusoid(radius, 0)
plt.plot(r, true, label='true')
plt.legend()









    Out[55]:





<matplotlib.legend.Legend at 0x11e42d3c8>



In [14]:

    
def linear(r, phi):
    return r

The same oscillatory behavior can be seen with a linear function



In [74]:

    
for index, zernike in enumerate(zernike_list):
    if index >= 6:
        zernike.generateCoefficients(linear, num_rintervals=10, num_aintervals=10,
                                 normalization="standard")
    

for r_index, radius in enumerate(r):
    for z_index, zernike in enumerate(zernike_list):
        if z_index >= 6:
            z_evals[z_index][r_index] = zernike(radius, 0, normalization="orthonormal")

for order, z in enumerate(z_evals):
    if order % 2 == 0 and order > 6:
        plt.plot(r, z, label=str(order))
plt.plot(r, r, label='true')
plt.legend()









    Out[74]:





<matplotlib.legend.Legend at 0x11eed5f60>

And again can be removed by increasing the accuracy of the integration



In [ ]:

    
for index, zernike in enumerate(zernike_list):
    if index >= 6:
        zernike.generateCoefficients(linear, num_rintervals=200, num_aintervals=200,
                                 normalization="standard")
    

for r_index, radius in enumerate(r):
    for z_index, zernike in enumerate(zernike_list):
        if z_index >= 6:
            z_evals[z_index][r_index] = zernike(radius, 0, normalization="orthonormal")



In [71]:

    
for order, z in enumerate(z_evals):
    if order % 2 == 0 and order > 6:
        plt.plot(r, z, label=str(order))
plt.plot(r, r, label='true')
plt.legend()









    Out[71]:





<matplotlib.legend.Legend at 0x11ecc4d68>

Tool for checking Zernike polynomial generation



In [36]:

    
import sympy as sp

def coeffs(m, n):
    clist = []
    clist.append(2 * (n - 1) * n * (n - 2))
    clist.append((n - 1) * (-m**2 - n * (n - 2)))
    clist.append(-n * (n + m - 2) * (n - m - 2) / 2)
    return [item / ((n + m) * (n - m) * (n - 2)) * 2 for item in clist]

rho = sp.symbols('rho')

def R(m, n):
    if m == n:
        return rho**n
    elif n == m + 2:
        return ((m + 2) * rho**2 - (m + 1)) * rho**m
    else:
        coeff1, coeff2, coeff3 = coeffs(m, n)
        return (coeff1 * rho**2 + coeff2) * R(m, n-2) + coeff3 * R(m, n-4)



In [ ]:

    
from series.Zernike import Zernike

max_order = 10
zernike_list = [Zernike(order, 0, 0, 60) for order in range(max_order + 1)]

zernike_tally.mean.shape = (zernike_tally.mean.size)

for order in range(max_order + 1):
    n = (order + 1) * (order + 2) // 2
    zernike_list[order].coefficients = zernike_tally.mean[:n]

import numpy as np

r = np.linspace(0, 60, 100)
z_evals = [np.zeros((r.size)) for zernike in zernike_list]

for r_index, radius in enumerate(r):
    for z_index, zernike in enumerate(zernike_list):
        z_evals[z_index][r_index] = zernike(radius, 0, normalization="orthonormal")



In [13]:

    
for order, z in enumerate(z_evals):
    if order % 2 == 0 and order <= 4:
        plt.plot(r, z, label=str(order))
plt.legend()









    Out[13]:





<matplotlib.legend.Legend at 0x1262d1278>

Using the code above, we determined that the radial Zernike polynomials in the MOOSE FET module are accurate!