<jupyter_start><jupyter_text>Lab 6 - Analysis of Heterogeneous ReactorsThis is a big lab. For the first time we venture outside the territory for which diffusion theory analysis can be used.  This lab includes, what I think to be, a good variety of pratical "experimentation" to investigate the impact of heterogenaity on reactor core physics.This was also an opportunity for students to use the MCNP plotting features (which I had been avoiding due to the mild complexity of using the MCNP plot system in a Windows environment).  The idea is to consider a reactor with light water moderator/coolant and low enriched uranium (metallic) fuel.  We will consider the following cases:a) homogeneous reactor with water, U-238 and U-235 all mixed together;b) A uranium metal fuel pin in a water tank (of the same size as the case above); the fuel pin is sized such that the *total* nuber of atoms of each isotope is unchanged; just arranged differently.c) An arrangement of 4 uranium metallic fuel pins; again the total number of atoms of each isotope is constant from the previous case;d) A 31x31 lattice of 2cm diameter pins (this will work out to be the same fuel mass and volume; this case is re-run with different lattice pitches. (2.5 cm; 3cm; 3.2cm; and 4cm - though I can pick more and more appropriate values for this analysis...)  Homogeneous mixture<jupyter_code>import openmc
import openmc.model
import numpy as np
import os<jupyter_output><empty_output><jupyter_text>Let's make a separte directory for each case so we can go back and look at the data<jupyter_code>if (not os.path.isdir('./homogeneous')): #
    os.mkdir('./homogeneous');
os.chdir('./homogeneous');<jupyter_output><empty_output><jupyter_text>Materials<jupyter_code>h_core = openmc.Material(name='core');
h_core.add_nuclide('U238',0.70340,'wo');
h_core.add_nuclide('U235',0.007105,'wo');
h_core.add_nuclide('H1',0.03216,'wo');
h_core.add_nuclide('O16',0.2573,'wo');
h_core.add_s_alpha_beta('c_H_in_H2O');
h_core.set_density('g/cm3',3.06);

mf = openmc.Materials([h_core]);
mf.export_to_xml();<jupyter_output><empty_output><jupyter_text>GeometryFor the homogeneous case, the geometry is just a cylindrical tank 155 cm in height and 75 cm in radius.<jupyter_code>core_tank = openmc.model.RightCircularCylinder([0.,0.,-5.],
                                              150.,radius=75.,
                                               axis='z',
                                              boundary_type='vacuum');

core_cell = openmc.Cell();
core_cell.fill = h_core;
core_cell.region = -core_tank;

root_universe = openmc.Universe();
root_universe.add_cells([core_cell]);

geometry = openmc.Geometry(root_universe);
geometry.export_to_xml();<jupyter_output><empty_output><jupyter_text>TalliesLet's set up some tallies to keep track of absorption/capture/etc... in the fuel and moderator.<jupyter_code>energy_filter = openmc.EnergyFilter([0.,0.625,20.e6]);
cell_filter = openmc.CellFilter([core_cell]);

tallies = openmc.Tallies()

t_abs = openmc.Tally(name='abs. tally');
t_abs.filters = [energy_filter,cell_filter];
t_abs.scores = ['absorption'];
t_abs.nuclides = ['H1','O16','U235','U238'];
tallies.append(t_abs);


t_fiss = openmc.Tally(name='fission tally');
t_fiss.filters = [energy_filter,cell_filter];
t_fiss.nuclides = ['U235','U238'];
t_fiss.scores = ['fission','nu-fission'];
tallies.append(t_fiss);

tallies.export_to_xml();<jupyter_output><empty_output><jupyter_text>Settings(will want to add tallies later...)<jupyter_code>settings = openmc.Settings();
settings.batches = 250;
settings.inactive = 100;
settings.particles = 10000;

bounds = [-75.,-75.,-5.,75.,75.,145.];
uniform_dist = openmc.stats.Box(bounds[:3],bounds[3:],
                               only_fissionable=True);
settings.source = openmc.source.Source(space=uniform_dist);

settings.export_to_xml();

openmc.run()<jupyter_output>%%%%%%%%%%%%%%%
                           %%%%%%%%%%%%%%%%%%%%%%%%
                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    %%%%%%%%%%%%%%%%%%%%%%%%
                                     %%%%%%%%%%%%%%%%%%%%%%%%
                 ###############      %%%%%%%%%%%%%%%%%%%%%%%%
                ##################     %%%%%%%%%%%%%%%%%%%%%%%
                ###################     %%%%%%%%%%%%%%%%%%%%%%%
                ####################     %%%%%%%%%%%%%%%%%%%%%%
                #####################     %%%%%%%%%%%%%%%%%%%%%
                ######################     %%%%%%%%%%%%%%%%%%%%
                #######################     %%%%%%%%%%%%%%%%%%
                 #######################     %%%%%%%%%%%%%%%%%
                 #####################[...]<jupyter_text>The reported $k_{\text{eff}}$ of 0.89210 is slightly different than my MCNP result (0.88996).Further investigation revealed that in the MCNP output to which I'm comparing, I didn't use the thermal scattering treatment for the light water.  If I take out the s_alpha_beta from the OpenMC model the revised k_eff is 0.89051 which is closer to the same result.  Lesson learned (again!): that shit matters!<jupyter_code>sp = openmc.StatePoint('statepoint.250.h5');
sp.tallies
abs_rate = sp.get_tally(name='abs. tally');
ar_df = abs_rate.get_pandas_dataframe();
ar_df.head(10)


total_abs = ar_df['mean'].sum();
print(f'Total absorption per source particle: %5.4f'%total_abs);
print(f'Leakage per source particle: %5.4f'%(1-total_abs));<jupyter_output>Total absorption per source particle: 0.9647
Leakage per source particle: 0.0353<jupyter_text>This leakage percent is close to what MCNP gives (3.91%) - but, remember, the MCNP input (for some reason) did not include thermal scattering treatment.  (if I re-try this without thermal scattering treatment: I get leakage per source particle of 3.85% which is closer to the MCNP result - consistent.)<jupyter_code>fiss_rate = sp.get_tally(name='fission tally');
fr_df = fiss_rate.get_pandas_dataframe();
fr_df.head(10)
fiss = fr_df[fr_df['score']=='fission'];
total_fission = fiss['mean'].sum()
print(f'Total fission: %5.4f per source particle'%total_fission);
fiss.head()<jupyter_output><empty_output><jupyter_text>I'm sure a lot more can be done here, but I want to move on to the next case, for now.<jupyter_code>os.chdir('..')
if (not os.path.isdir('./one_pin')): #
    os.mkdir('./one_pin');
os.chdir('./one_pin');<jupyter_output><empty_output><jupyter_text>One PinHere I will collect all of the uranium in the homogeneous core mixture and combine into a single pin.  This pin will be placed into a tank of water the same size as before  Material<jupyter_code>fuel = openmc.Material(name='fuel');
fuel.add_nuclide('U238',0.99,'wo');
fuel.add_nuclide('U235',0.01,'wo');
fuel.set_density('g/cm3',19.1);

mod = openmc.Material(name='mod');
mod.add_nuclide('H1',2,'ao');
mod.add_nuclide('O16',1,'ao');
mod.add_s_alpha_beta('c_H_in_H2O');
mod.set_density('g/cm3',1.0);

mf = openmc.Materials();
mf += [fuel,mod];

mf.export_to_xml();<jupyter_output><empty_output><jupyter_text>GeometryThe fuel pin is a rcc 100 cm high and 31 cm in radius.  This sits inside the larger cylindrical tank (same dimensions of the homogeneous core) already defined as "core_tank"<jupyter_code>fuel_pin = openmc.model.RightCircularCylinder([0.,0.,0.],
                                              100.,radius=31.,
                                               axis='z');

fuel_cell = openmc.Cell();
fuel_cell.fill = fuel;
fuel_cell.region = -fuel_pin;

rx_tank = openmc.Cell();
rx_tank.fill = mod;
rx_tank.region = +fuel_pin & -core_tank;

root = openmc.Universe();
root.add_cells([fuel_cell,rx_tank]);

geometry = openmc.Geometry();
geometry.root_universe = root;

geometry.export_to_xml();<jupyter_output><empty_output><jupyter_text>TalliesI will use the same tallies but I'll need to change the cell filters for the current geometry so some adjustment to the objects will be required.<jupyter_code>#energy_filter = openmc.EnergyFilter([0.,0.625,20.e6]); #<-- keep this for reference

fuel_cell_filter = openmc.CellFilter([fuel_cell]);
rx_tank_cell_filter = openmc.CellFilter([rx_tank]);

tallies = openmc.Tallies()

abs_fuel = openmc.Tally(name='abs. fuel');
abs_fuel.filters = [energy_filter,fuel_cell_filter];
abs_fuel.scores = ['absorption'];
abs_fuel.nuclides = ['U235','U238'];
tallies.append(abs_fuel);

abs_mod = openmc.Tally(name='abs. mod');
abs_mod.filters = [energy_filter,rx_tank_cell_filter];
abs_mod.nuclides = ['H1','O16'];
abs_mod.scores = ['absorption'];
tallies.append(abs_mod);



t_fiss = openmc.Tally(name='fission tally');
t_fiss.filters = [energy_filter,fuel_cell_filter];
t_fiss.nuclides = ['U235','U238'];
t_fiss.scores = ['fission','nu-fission'];
tallies.append(t_fiss);

tallies.export_to_xml();<jupyter_output><empty_output><jupyter_text>SettingsI will use mostly the same settings.  I just want to change the bounding box so it now only encompases the fuel region.  Probably not strictly required but anyway...<jupyter_code>bounds = [-31.,-31.,0.,31.,31.,100.];
uniform_dist = openmc.stats.Box(bounds[:3],bounds[3:],
                               only_fissionable=True);
settings.source = openmc.source.Source(space=uniform_dist);

settings.export_to_xml();
openmc.run()<jupyter_output>%%%%%%%%%%%%%%%
                           %%%%%%%%%%%%%%%%%%%%%%%%
                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    %%%%%%%%%%%%%%%%%%%%%%%%
                                     %%%%%%%%%%%%%%%%%%%%%%%%
                 ###############      %%%%%%%%%%%%%%%%%%%%%%%%
                ##################     %%%%%%%%%%%%%%%%%%%%%%%
                ###################     %%%%%%%%%%%%%%%%%%%%%%%
                ####################     %%%%%%%%%%%%%%%%%%%%%%
                #####################     %%%%%%%%%%%%%%%%%%%%%
                ######################     %%%%%%%%%%%%%%%%%%%%
                #######################     %%%%%%%%%%%%%%%%%%
                 #######################     %%%%%%%%%%%%%%%%%
                 #####################[...]<jupyter_text>The result is much closer to the MCNP result (why?).  The main difference from the homgeneous case is an increase in neutron capture.  SThe tallies will now be analyzed to determine:a) where does most of the neutron capture occur;b) in which nuclide; andc) at what energy<jupyter_code>sp = openmc.StatePoint('statepoint.250.h5');
sp.tallies
abs_fuel = sp.get_tally(name='abs. fuel');
abs_fuel_df = abs_fuel.get_pandas_dataframe();
abs_fuel_df.head(10)
abs_mod = sp.get_tally(name='abs. mod');
abs_mod_df = abs_mod.get_pandas_dataframe();
abs_mod_df.head(10)<jupyter_output><empty_output><jupyter_text>A brief glance at the above results indicate that the vast majority of the absorptions are happeinging:a) epithermal and fast neutrons in U238; andb) thermal neutrons in H1 in the moderator.<jupyter_code>os.chdir('..')
if (not os.path.isdir('./four_pins')): #
    os.mkdir('./four_pins');
os.chdir('./four_pins');<jupyter_output><empty_output><jupyter_text>4-PinsIn this iteration we divide the 1 large pin into 4 equally sized pins that, cumulatively, has the same volume as the single pin.  In making this change we should see a reduction in capture in U-238 but - well, let's just do it and find out for sure.  While we're at it, we will practice making lattices in OpenMC, practice plotting the geometry, and learn how to do tallies in lattices.  MaterialsThere are no changes to the materials for this part of the project so we will merely export the existing materials object to xml.<jupyter_code># no changes; create a local copy of the xml file.
mf.export_to_xml();

# reminder the two materials are:
# fuel - 1% enriched uranium fuel 
# mod - light water moderator<jupyter_output><empty_output><jupyter_text>GeometryI will make a 2x2 lattice of the fuel pins.  The lattice pitch is 38 cm; the entire area is in a RPP centered on the origin that is 76 cm square and 100 cm tall.<jupyter_code>pitch = 38; #cm
fuel_radius = 15.5; #cm

# create a RCC for a generic fuel pin in the fuel pin universe

fuel_pin = openmc.model.RightCircularCylinder([0.,0.,0.],100.,
                                              radius=fuel_radius,
                                              axis='z');

fp_fuel = openmc.Cell(name='fuel');
fp_fuel.fill = fuel;
fp_fuel.region = -fuel_pin;

fp_mod = openmc.Cell(name='mod in fuel');
fp_mod.fill = mod;
fp_mod.region = +fuel_pin;

fp = openmc.Universe();
fp.add_cells((fp_fuel,fp_mod));

mod_cell = openmc.Cell();
mod_cell.fill = mod;

all_water = openmc.Universe();
all_water.add_cell(mod_cell);



colors = {}
colors[mod] = 'blue';
colors[fuel] = 'yellow';
# let's make some other plots
main = openmc.Cell();
main.fill = fp;

root = openmc.Universe();
root.add_cell(main);

g = openmc.Geometry();
g.root_universe = root;
g.export_to_xml();

p = openmc.Plot();
p.width = [pitch,pitch];
p.pixels = [400,400];
p.origin = [0.,0.,10.];
#p.basis = 'xy';
p.color_by = 'material';
p.colors = colors

openmc.plot_inline(p);
# maybe you don't really need this?
core_area = openmc.model.RectangularParallelepiped(-38.,38.,-38.,38.,
                                                  0.,100.);

lattice = openmc.RectLattice();
lattice.dimension = [2,2];
lattice.lower_left = [-pitch,-pitch];
lattice.pitch = [pitch,pitch];

lattice.universes = [
    [fp,fp],
    [fp,fp]
];

lattice.outer = all_water;

core_cell = openmc.Cell();
core_cell.fill = lattice;
core_cell.region = -core_area;

refl_cell = openmc.Cell();
refl_cell.fill = mod;
refl_cell.region = +core_area & -core_tank;



#main = openmc.Cell()
#main.fill = lattice;

root = openmc.Universe();
root.add_cells([core_cell,refl_cell]);

g = openmc.Geometry();
g.root_universe = root;
g.export_to_xml();<jupyter_output><empty_output><jupyter_text>I think this should all "work" -- even if it's not the best thing ever written.  To be sure, it's a good idea to make a plot to be sure<jupyter_code>colors = {}
colors[mod] = 'blue';
colors[fuel] = 'yellow';
p = openmc.Plot();
p.width = [4*pitch,4*pitch];
p.origin = [0.,0.,10.];
p.pixels = [400,400];
p.color_by = 'material';
p.colors = colors

openmc.plot_inline(p);<jupyter_output><empty_output><jupyter_text>Okay...this looks correct now.<jupyter_code>fuel_cell_filter = openmc.DistribcellFilter([fp_fuel]);
mod_in_fuel_filter = openmc.DistribcellFilter([fp_mod]);
rfl_cell_filter = openmc.CellFilter([refl_cell]);

tallies = openmc.Tallies()

abs_fuel = openmc.Tally(name='abs. fuel');
abs_fuel.filters = [energy_filter,fuel_cell_filter];
abs_fuel.scores = ['absorption'];
abs_fuel.nuclides = ['U235','U238'];
tallies.append(abs_fuel);

abs_mod_ref = openmc.Tally(name='abs. mod refl');
abs_mod_ref.filters = [energy_filter,rfl_cell_filter];
abs_mod_ref.nuclides = ['H1','O16'];
abs_mod_ref.scores = ['absorption'];
tallies.append(abs_mod_ref);

abs_mod_fuel = openmc.Tally(name='abs. mod fuel');
abs_mod_fuel.filters = [energy_filter,mod_in_fuel_filter];
abs_mod_fuel.nuclides = ['H1','O16'];
abs_mod_fuel.scores = ['absorption'];
tallies.append(abs_mod_fuel);


t_fiss = openmc.Tally(name='fission tally');
t_fiss.filters = [energy_filter,fuel_cell_filter];
t_fiss.nuclides = ['U235','U238'];
t_fiss.scores = ['fission','nu-fission'];
tallies.append(t_fiss);

tallies.export_to_xml();

settings.export_to_xml();

openmc.run();<jupyter_output>%%%%%%%%%%%%%%%
                           %%%%%%%%%%%%%%%%%%%%%%%%
                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    %%%%%%%%%%%%%%%%%%%%%%%%
                                     %%%%%%%%%%%%%%%%%%%%%%%%
                 ###############      %%%%%%%%%%%%%%%%%%%%%%%%
                ##################     %%%%%%%%%%%%%%%%%%%%%%%
                ###################     %%%%%%%%%%%%%%%%%%%%%%%
                ####################     %%%%%%%%%%%%%%%%%%%%%%
                #####################     %%%%%%%%%%%%%%%%%%%%%
                ######################     %%%%%%%%%%%%%%%%%%%%
                #######################     %%%%%%%%%%%%%%%%%%
                 #######################     %%%%%%%%%%%%%%%%%
                 #####################[...]<jupyter_text>Once again; this is much like the MCNP result.<jupyter_code>sp = openmc.StatePoint('statepoint.250.h5');
sp.tallies

abs_fuel = sp.get_tally(name='abs. fuel');
abs_fuel_df = abs_fuel.get_pandas_dataframe();
#abs_fuel_df.head(25)
u235_abs = abs_fuel_df[(abs_fuel_df['nuclide']=='U235')]['mean'].sum()
print(f'Absorptioni in U235: %5.4f'%u235_abs);
u238_abs = abs_fuel_df[abs_fuel_df['nuclide']=='U238']['mean'].sum()
print(f'Absorptioni in U238: %5.4f'%u238_abs);
abs_mod_refl = sp.get_tally(name='abs. mod refl');
abs_mod_refl_df = abs_mod_refl.get_pandas_dataframe();
abs_mod_refl_df.head(25)
total_abs_mod_refl = abs_mod_refl_df['mean'].sum();
print(f'total absorption in moderator in reflector region: %5.4f'%total_abs_mod_refl);
abs_mod_fuel = sp.get_tally(name='abs. mod fuel');
abs_mod_fuel_df = abs_mod_fuel.get_pandas_dataframe();
#abs_mod_fuel_df.head(25)
total_abs_mod_fuel = abs_mod_fuel_df['mean'].sum()
print(f'total absorption in the moderator in the core region: %5.4f'%total_abs_mod_fuel);
total_mod_abs = total_abs_mod_fuel+total_abs_mod_refl;
print(f'total absorption in moderator material: %5.4f'%total_mod_abs);
fiss = sp.get_tally(name='fission tally');
fiss_df = fiss.get_pandas_dataframe();
#fiss_df.head(25)
fission = fiss_df[fiss_df['score']=='fission'];
#fission.head(10)
fiss_235 = fission[fission['nuclide']=='U235'];
fiss_238 = fission[fission['nuclide']=='U238'];
fiss_235.head(10)

fiss_238.head(10)

print(f'fission in U238: %5.4f'%fiss_238['mean'].sum());
print(f'fission in U235: %5.4f'%fiss_235['mean'].sum());
os.chdir('..')
if (not os.path.isdir('./31x31_pins')): #
    os.mkdir('./31x31_pins');
os.chdir('./31x31_pins');<jupyter_output><empty_output><jupyter_text>31x31 latticesFor this section all of the materials, tallies, and settings can remain the same; the geometry will need to be updated to conver the 2x2 lattice of thick rods into a 31x31 lattice of much thinner rods of the same total volume.  We will start by defining a function that creates a new geometry with user input for pitch<jupyter_code>def gen_lattice(l_pitch):
    #input lattice pitch
    
    n_pins = 31;
    pin_r = 1.; # cm, fuel pin radius
    
    fuel_pin = openmc.model.RightCircularCylinder([0.,0.,0.],100.,
                                              radius=pin_r,
                                              axis='z');

    fp_fuel = openmc.Cell(name='fuel');
    fp_fuel.fill = fuel;
    fp_fuel.region = -fuel_pin;

    fp_mod = openmc.Cell(name='mod in fuel');
    fp_mod.fill = mod;
    fp_mod.region = +fuel_pin;

    fp = openmc.Universe();
    fp.add_cells((fp_fuel,fp_mod));

    mod_cell = openmc.Cell();
    mod_cell.fill = mod;

    all_water = openmc.Universe();
    all_water.add_cell(mod_cell);
    core_bd = l_pitch*n_pins/2.
    core_area = openmc.model.RectangularParallelepiped(-core_bd,core_bd,-core_bd,core_bd,
                                                  0.,100.);

    lattice = openmc.RectLattice();
    lattice.dimension = [n_pins,n_pins];
    lattice.lower_left = [-core_bd,-core_bd];
    lattice.pitch = [l_pitch,l_pitch];

    lattice.universes = [[ fp for i in range(n_pins)]for j in range(n_pins)];
           

    lattice.outer = all_water;

    core_cell = openmc.Cell();
    core_cell.fill = lattice;
    core_cell.region = -core_area;

    refl_cell = openmc.Cell();
    refl_cell.fill = mod;
    refl_cell.region = +core_area & -core_tank;

    root = openmc.Universe();
    root.add_cells([core_cell,refl_cell]);

    g = openmc.Geometry();
    g.root_universe = root;
    g.export_to_xml();
    
    
    fuel_cell_filter = openmc.DistribcellFilter([fp_fuel]);
    mod_in_fuel_filter = openmc.DistribcellFilter([fp_mod]);
    rfl_cell_filter = openmc.CellFilter([refl_cell]);

    tallies = openmc.Tallies()

    abs_fuel = openmc.Tally(name='abs. fuel');
    abs_fuel.filters = [energy_filter,fuel_cell_filter];
    abs_fuel.scores = ['absorption'];
    abs_fuel.nuclides = ['U235','U238'];
    tallies.append(abs_fuel);

    abs_mod_ref = openmc.Tally(name='abs. mod refl');
    abs_mod_ref.filters = [energy_filter,rfl_cell_filter];
    abs_mod_ref.nuclides = ['H1','O16'];
    abs_mod_ref.scores = ['absorption'];
    tallies.append(abs_mod_ref);

    abs_mod_fuel = openmc.Tally(name='abs. mod fuel');
    abs_mod_fuel.filters = [energy_filter,mod_in_fuel_filter];
    abs_mod_fuel.nuclides = ['H1','O16'];
    abs_mod_fuel.scores = ['absorption'];
    tallies.append(abs_mod_fuel);


    t_fiss = openmc.Tally(name='fission tally');
    t_fiss.filters = [energy_filter,fuel_cell_filter];
    t_fiss.nuclides = ['U235','U238'];
    t_fiss.scores = ['fission','nu-fission'];
    tallies.append(t_fiss);

    tallies.export_to_xml();
    
# generate the first geometry and make a plot
if (not os.path.isdir('./p2_5')): #
    os.mkdir('./p2_5');
os.chdir('./p2_5');


mf.export_to_xml();
settings.export_to_xml();
gen_lattice(2.5);
# do without extra tallies for now.

p = openmc.Plot();
p.width = [200., 200.];
p.origin = [0.,0.,10.];
p.pixels = [1024,1024];
p.color_by = 'material';
p.colors = colors

openmc.plot_inline(p);
openmc.run();
os.chdir('..')
if (not os.path.isdir('./p3_0')): #
    os.mkdir('./p3_0');
os.chdir('./p3_0');


mf.export_to_xml();
settings.export_to_xml();
gen_lattice(3.0);
# do without extra tallies for now.

p = openmc.Plot();
p.width = [200., 200.];
p.origin = [0.,0.,10.];
p.pixels = [1024,1024];
p.color_by = 'material';
p.colors = colors

openmc.plot_inline(p);
openmc.run()
os.chdir('..')
if (not os.path.isdir('./p3_2')): #
    os.mkdir('./p3_2');
os.chdir('./p3_2');


mf.export_to_xml();
settings.export_to_xml();
gen_lattice(3.2);
# do without extra tallies for now.

p = openmc.Plot();
p.width = [200., 200.];
p.origin = [0.,0.,10.];
p.pixels = [1024,1024];
p.color_by = 'material';
p.colors = colors

openmc.plot_inline(p);
openmc.run();
os.chdir('..')
if (not os.path.isdir('./p3_4')): #
    os.mkdir('./p3_4');
os.chdir('./p3_4');


mf.export_to_xml();
settings.export_to_xml();
gen_lattice(3.4);
# do without extra tallies for now.

p = openmc.Plot();
p.width = [200., 200.];
p.origin = [0.,0.,10.];
p.pixels = [1024,1024];
p.color_by = 'material';
p.colors = colors

openmc.plot_inline(p);
openmc.run()
os.chdir('../..')
if (not os.path.isdir('./partB')): #
    os.mkdir('./partB');
os.chdir('./partB');<jupyter_output><empty_output><jupyter_text>Part BThis will be the same as the 31x31 lattice except the core region will use heavy water and the reflector region will be graphite.  Materials<jupyter_code>fuel = openmc.Material(name='fuel');
fuel.add_nuclide('U238',0.99,'wo');
fuel.add_nuclide('U235',0.01,'wo');
fuel.set_density('g/cm3',19.1);

mod = openmc.Material(name='mod');
mod.add_nuclide('H2',2,'ao');
mod.add_nuclide('O16',1,'ao');
mod.add_s_alpha_beta('c_D_in_D2O');
mod.set_density('g/cm3',1.105);

graph_refl = openmc.Material(name='graph_refl');
graph_refl.add_nuclide('C0',1,'ao');
graph_refl.add_s_alpha_beta('c_Graphite');
graph_refl.set_density('g/cm3',1.7);



mf = openmc.Materials();
mf += [fuel,mod,graph_refl];

mf.export_to_xml();<jupyter_output><empty_output><jupyter_text>Geometry<jupyter_code>def gen_lattice_hwg(l_pitch):
    #input lattice pitch
    
    n_pins = 31;
    pin_r = 1.; # cm, fuel pin radius
    
    fuel_pin = openmc.model.RightCircularCylinder([0.,0.,0.],100.,
                                              radius=pin_r,
                                              axis='z');

    fp_fuel = openmc.Cell(name='fuel');
    fp_fuel.fill = fuel;
    fp_fuel.region = -fuel_pin;

    fp_mod = openmc.Cell(name='mod in fuel');
    fp_mod.fill = mod;
    fp_mod.region = +fuel_pin;

    fp = openmc.Universe();
    fp.add_cells((fp_fuel,fp_mod));

    mod_cell = openmc.Cell();
    mod_cell.fill = mod;

    all_water = openmc.Universe();
    all_water.add_cell(mod_cell);
    
    core_bd = l_pitch*n_pins/2.
    core_area = openmc.model.RectangularParallelepiped(-core_bd,core_bd,-core_bd,core_bd,
                                                  0.,100.);

    lattice = openmc.RectLattice();
    lattice.dimension = [n_pins,n_pins];
    lattice.lower_left = [-core_bd,-core_bd];
    lattice.pitch = [l_pitch,l_pitch];

    lattice.universes = [[ fp for i in range(n_pins)]for j in range(n_pins)];
           

    lattice.outer = all_water;

    core_cell = openmc.Cell();
    core_cell.fill = lattice;
    core_cell.region = -core_area;

    refl_cell = openmc.Cell();
    refl_cell.fill = graph_refl;
    refl_cell.region = +core_area & -core_tank;

    root = openmc.Universe();
    root.add_cells([core_cell,refl_cell]);

    g = openmc.Geometry();
    g.root_universe = root;
    g.export_to_xml();
    
    
    fuel_cell_filter = openmc.DistribcellFilter([fp_fuel]);
    mod_in_fuel_filter = openmc.DistribcellFilter([fp_mod]);
    rfl_cell_filter = openmc.CellFilter([refl_cell]);

    tallies = openmc.Tallies()

    abs_fuel = openmc.Tally(name='abs. fuel');
    abs_fuel.filters = [energy_filter,fuel_cell_filter];
    abs_fuel.scores = ['absorption'];
    abs_fuel.nuclides = ['U235','U238'];
    tallies.append(abs_fuel);

    abs_mod_ref = openmc.Tally(name='abs. mod refl');
    abs_mod_ref.filters = [energy_filter,rfl_cell_filter];
    abs_mod_ref.nuclides = ['C0'];
    abs_mod_ref.scores = ['absorption'];
    tallies.append(abs_mod_ref);

    abs_mod_fuel = openmc.Tally(name='abs. mod fuel');
    abs_mod_fuel.filters = [energy_filter,mod_in_fuel_filter];
    abs_mod_fuel.nuclides = ['H2','O16'];
    abs_mod_fuel.scores = ['absorption'];
    tallies.append(abs_mod_fuel);


    t_fiss = openmc.Tally(name='fission tally');
    t_fiss.filters = [energy_filter,fuel_cell_filter];
    t_fiss.nuclides = ['U235','U238'];
    t_fiss.scores = ['fission','nu-fission'];
    tallies.append(t_fiss);

    tallies.export_to_xml();
gen_lattice_hwg(3.2); #lattice with heavy water and graphite, 3.2 cm pitch
settings.export_to_xml();
colors = {}
colors[mod] = 'blue';
colors[fuel] = 'yellow';
colors[graph_refl] = 'gray';
p = openmc.Plot();
p.width = [200., 200.];
p.origin = [0.,0.,10.];
p.pixels = [1024,1024];
p.color_by = 'material';
p.colors = colors

openmc.plot_inline(p);
openmc.run()
os.chdir('..')
if (not os.path.isdir('./partC')): #
    os.mkdir('./partC');
os.chdir('./partC');<jupyter_output><empty_output><jupyter_text>Part CIn this iteration of the analysis we increase the overall size of the reactor tank to 300 cm and, using the same 31x31 lattice, vary the pitch from 4cm to 12cm noting the change of $k_{\text{eff}}$ along the way.    Materials, Geometry, and TalliesSo - a few small changes are required but I thought it would be nice to have all of the relevent bits of code all together (for once).  T<jupyter_code>fuel = openmc.Material(name='fuel');
fuel.add_nuclide('U238',0.99,'wo');
fuel.add_nuclide('U235',0.01,'wo');
fuel.set_density('g/cm3',19.1);

mod = openmc.Material(name='mod');
mod.add_nuclide('H2',2,'ao');
mod.add_nuclide('O16',1,'ao');
mod.add_s_alpha_beta('c_D_in_D2O');
mod.set_density('g/cm3',1.105);



mf = openmc.Materials();
mf += [fuel,mod];

mf.export_to_xml();


core_tank = openmc.model.RightCircularCylinder([0.,0.,-5.],
                                              150.,radius=300.,
                                               axis='z',
                                              boundary_type='vacuum');

core_cell = openmc.Cell();
core_cell.fill = h_core;
core_cell.region = -core_tank;

root_universe = openmc.Universe();
root_universe.add_cells([core_cell]);

geometry = openmc.Geometry(root_universe);
geometry.export_to_xml();


def gen_lattice_hw(l_pitch):
    #input lattice pitch
    
    n_pins = 31;
    pin_r = 1.; # cm, fuel pin radius
    
    fuel_pin = openmc.model.RightCircularCylinder([0.,0.,0.],100.,
                                              radius=pin_r,
                                              axis='z');

    fp_fuel = openmc.Cell(name='fuel');
    fp_fuel.fill = fuel;
    fp_fuel.region = -fuel_pin;

    fp_mod = openmc.Cell(name='mod in fuel');
    fp_mod.fill = mod;
    fp_mod.region = +fuel_pin;

    fp = openmc.Universe();
    fp.add_cells((fp_fuel,fp_mod));

    mod_cell = openmc.Cell();
    mod_cell.fill = mod;

    all_water = openmc.Universe();
    all_water.add_cell(mod_cell);
    
    core_bd = l_pitch*n_pins/2.
    core_area = openmc.model.RectangularParallelepiped(-core_bd,core_bd,-core_bd,core_bd,
                                                  0.,100.);

    lattice = openmc.RectLattice();
    lattice.dimension = [n_pins,n_pins];
    lattice.lower_left = [-core_bd,-core_bd];
    lattice.pitch = [l_pitch,l_pitch];

    lattice.universes = [[ fp for i in range(n_pins)]for j in range(n_pins)];
           

    lattice.outer = all_water;

    core_cell = openmc.Cell();
    core_cell.fill = lattice;
    core_cell.region = -core_area;

    refl_cell = openmc.Cell();
    refl_cell.fill = mod;
    refl_cell.region = +core_area & -core_tank;

    root = openmc.Universe();
    root.add_cells([core_cell,refl_cell]);

    g = openmc.Geometry();
    g.root_universe = root;
    g.export_to_xml();
    
    
    fuel_cell_filter = openmc.DistribcellFilter([fp_fuel]);
    mod_in_fuel_filter = openmc.DistribcellFilter([fp_mod]);
    rfl_cell_filter = openmc.CellFilter([refl_cell]);

    tallies = openmc.Tallies()

    abs_fuel = openmc.Tally(name='abs. fuel');
    abs_fuel.filters = [energy_filter,fuel_cell_filter];
    abs_fuel.scores = ['absorption'];
    abs_fuel.nuclides = ['U235','U238'];
    tallies.append(abs_fuel);

    abs_mod_ref = openmc.Tally(name='abs. mod refl');
    abs_mod_ref.filters = [energy_filter,rfl_cell_filter];
    abs_mod_ref.nuclides = ['H2','O16'];
    abs_mod_ref.scores = ['absorption'];
    tallies.append(abs_mod_ref);

    abs_mod_fuel = openmc.Tally(name='abs. mod fuel');
    abs_mod_fuel.filters = [energy_filter,mod_in_fuel_filter];
    abs_mod_fuel.nuclides = ['H2','O16'];
    abs_mod_fuel.scores = ['absorption'];
    tallies.append(abs_mod_fuel);


    t_fiss = openmc.Tally(name='fission tally');
    t_fiss.filters = [energy_filter,fuel_cell_filter];
    t_fiss.nuclides = ['U235','U238'];
    t_fiss.scores = ['fission','nu-fission'];
    tallies.append(t_fiss);

    tallies.export_to_xml();
gen_lattice_hw(4.0);

p = openmc.Plot();
p.width = [400., 400.];
p.origin = [0.,0.,10.];
p.pixels = [1024,1024];
p.color_by = 'material';
p.colors = colors

openmc.plot_inline(p);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.188538553016424, 'fraction_numerical': 0.024587792884003472, 'mean_word_length': 3.6335070737155624, 'pattern_counts': {'":': 0, '<': 94, '<?xml version=': 0, '>': 93, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 33}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 46, 'profanity_word_list.json': 1, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '24', 'dir': 'jupyter-structured-clean-dedup', 'id': '190550', 'license': 'MIT', 'n_tokens_mistral': 11522, 'n_tokens_neox': 10485, 'n_words': 2132, 'path': 'lab6/.ipynb_checkpoints/lab6-checkpoint.ipynb', 'repo_name': 'stu314159/er362_omc', 'size': '31113'}