<jupyter_start><jupyter_text>Part 2 - Simulations of effective dose as a function of distanceSometimes it is useful to have the dose as function of distance from the source. The following example builds on the previous example by finding the dose as a function of distance.The geometry is a little more complex in this example as we need additional surfaces to tally on. The following section makes the geometry, assigns materials and plots the result.This is the similar to the previous task geometry, but this time we have nested spheres which will be used to tally the dose on.<jupyter_code>import openmc

# dimentions are in cm
height = 150
outer_radius = 74.5 / 2
thickness = 1.5

outer_cylinder = openmc.ZCylinder(r=outer_radius)
inner_cylinder = openmc.ZCylinder(r=outer_radius-thickness)
inner_top = openmc.ZPlane(z0=height*0.5)
inner_bottom = openmc.ZPlane(z0=-height*0.5)
outer_top = openmc.ZPlane(z0=(height*0.5)+thickness)
outer_bottom = openmc.ZPlane(z0=(-height*0.5)-thickness)

# these spheres are used to tally the dose
sphere_1 = openmc.Sphere(r=200)
sphere_2 = openmc.Sphere(r=250)
sphere_3 = openmc.Sphere(r=300)
sphere_4 = openmc.Sphere(r=350)
sphere_5 = openmc.Sphere(r=400)
sphere_6 = openmc.Sphere(r=450)
sphere_7 = openmc.Sphere(r=500)
# can't actually tally on the end of universe sphere so sphere_8 is not in the tally
sphere_8 = openmc.Sphere(r=501, boundary_type='vacuum')

steel = openmc.Material(name='steel')
steel.set_density('g/cm3', 7.75)
steel.add_element('Fe', 0.95, percent_type='wo')
steel.add_element('C', 0.05, percent_type='wo')

mats = openmc.Materials([steel])

cylinder_region = -outer_cylinder & +inner_cylinder & -inner_top & +inner_bottom
cylinder_cell = openmc.Cell(region=cylinder_region)
cylinder_cell.fill = steel

top_cap_region = -outer_top & +inner_top & -outer_cylinder
top_cap_cell = openmc.Cell(region=top_cap_region)
top_cap_cell.fill = steel

bottom_cap_region = +outer_bottom & -inner_bottom & -outer_cylinder
bottom_cap_cell = openmc.Cell(region=bottom_cap_region)
bottom_cap_cell.fill = steel

inner_void_region = -inner_cylinder & -inner_top & +inner_bottom
inner_void_cell = openmc.Cell(region=inner_void_region)

# sphere 1 region is below -sphere_1 and not (~) in the other regions
sphere_1_region = -sphere_1
sphere_1_cell = openmc.Cell(
    region= sphere_1_region
    & ~bottom_cap_region
    & ~top_cap_region
    & ~cylinder_region
    & ~inner_void_region
)

sphere_2_region = +sphere_1 & -sphere_2
sphere_2_cell = openmc.Cell(region= sphere_2_region)

sphere_3_region = +sphere_2 & -sphere_3
sphere_3_cell = openmc.Cell(region= sphere_3_region)

sphere_4_region = +sphere_3 & -sphere_4
sphere_4_cell = openmc.Cell(region= sphere_4_region)

sphere_5_region = +sphere_4 & -sphere_5
sphere_5_cell = openmc.Cell(region= sphere_5_region)

sphere_6_region = +sphere_5 & -sphere_6
sphere_6_cell = openmc.Cell(region= sphere_6_region)

sphere_7_region = +sphere_6 & -sphere_7
sphere_7_cell = openmc.Cell(region= sphere_7_region)

sphere_8_region = +sphere_7 & -sphere_8
sphere_8_cell = openmc.Cell(region= sphere_8_region)

universe = openmc.Universe(cells=[
    inner_void_cell, cylinder_cell, top_cap_cell,
    bottom_cap_cell, sphere_1_cell, sphere_2_cell,
    sphere_2_cell, sphere_3_cell, sphere_4_cell,
    sphere_5_cell, sphere_6_cell, sphere_7_cell,
    sphere_8_cell])

geom = openmc.Geometry(universe)<jupyter_output><empty_output><jupyter_text>This section plots the geometry and colours it according to cell.<jupyter_code>import matplotlib.pyplot as plt

color_assignment = {sphere_1_cell: 'green',
                    sphere_2_cell: 'brown',
                    sphere_3_cell: 'lime',
                    sphere_4_cell: 'lightgray',
                    sphere_5_cell: 'maroon',
                    sphere_6_cell: 'magenta',
                    sphere_7_cell: 'cyan',
                    sphere_8_cell: 'purple',
                    inner_void_cell: 'grey',
                    bottom_cap_cell: 'red',
                    top_cap_cell: 'blue',
                    cylinder_cell:'yellow',
                   }
x, y = 1000, 1000
plt.show(universe.plot(width=(x, y), basis='xz', colors=color_assignment))
plt.show(universe.plot(width=(x, y), basis='xy', colors=color_assignment))<jupyter_output><empty_output><jupyter_text>This section makes the source, note the use of the Co60 gamma source with two energy levels.<jupyter_code># Instantiate a Settings object
sett = openmc.Settings()
sett.batches = 10
sett.inactive = 0
sett.particles = 500
sett.run_mode = 'fixed source'

# Create a gamma point source
source = openmc.Source()
source.space = openmc.stats.Point((0, 0, 0))
source.angle = openmc.stats.Isotropic()
# This is a Co60 source, see the task on sources to understand it
source.energy = openmc.stats.Discrete([1.1732e6,1.3325e6], [0.5, 0.5])
source.particle = 'photon'
sett.source = source

energy_bins_p, dose_coeffs_p = openmc.data.dose_coefficients(
    particle='photon',
    geometry='AP'
)

energy_function_filter_p = openmc.EnergyFunctionFilter(
    energy_bins_p,
    dose_coeffs_p
)

photon_particle_filter = openmc.ParticleFilter(["photon"])

tallies = openmc.Tallies()

surfaces_to_tally = [
    sphere_1, sphere_2, sphere_3, sphere_4, sphere_5, sphere_6, sphere_7
]

# this loops adds tallies for each sphere surface
for surface_id, surface in zip(range(1,9), surfaces_to_tally):
    surface_filter = openmc.SurfaceFilter(surface)
    dose_tally = openmc.Tally(name="dose_tally_on_surface_"+str(surface_id))
    dose_tally.scores = ["current"]
    dose_tally.filters = [
        surface_filter,
        photon_particle_filter,
        energy_function_filter_p,
    ]
    tallies.append(dose_tally)<jupyter_output><empty_output><jupyter_text>This section runs the simulation.<jupyter_code># Run OpenMC!
model = openmc.model.Model(geom, mats, sett, tallies)
!rm *.h5
sp_filename = model.run()<jupyter_output><empty_output><jupyter_text>This section loads up the results and post processes them so that the units are more useful.<jupyter_code>import math

# open the results file
sp = openmc.StatePoint(sp_filename)

dose_rates_in_pSv = []
distances = []
# access the tally using pandas dataframes
for surface_id, surface in zip(range(1,9), surfaces_to_tally):
    tally = sp.get_tally(name='dose_tally_on_surface_'+str(surface_id))
    df = tally.get_pandas_dataframe()
    tally_result = df['mean'].sum()
    tally_std_dev = df['std. dev.'].sum()

    # convert from the tally output units of pSv cm² to pSv by dividing by the surface area of the surface
    dose_in_pSv = tally_result / (4 * math.pi * math.pow(surface.r, 2))
    distances.append(surface.r)

    source_activity = 56000  # in decays per second (Bq)
    emission_rate = 2  # the number of gammas emitted per decay which is approximately 2 for Co60
    gamma_per_second = source_activity * emission_rate
    dose_rate_in_pSv = dose_in_pSv * gamma_per_second
    
    dose_rates_in_pSv.append(dose_rate_in_pSv)

# print results
print('The surface dose = ', dose_rates_in_pSv, 'pSv per second')
print('At distances = ', distances)<jupyter_output><empty_output><jupyter_text>This section plots the dose as a function of distance from the source.<jupyter_code>import matplotlib.pyplot as plt
plt.plot(distances, dose_rates_in_pSv, 'o')
plt.ylabel('Dose rates (pSv per second)')
plt.xlabel('Distance (cm)')
plt.savefig('dose_rate_from_co60_cask.svg')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09884346422807962, 'fraction_numerical': 0.024878967186659495, 'mean_word_length': 4.672768878718536, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '145740', 'license': 'MIT', 'n_tokens_mistral': 2762, 'n_tokens_neox': 2582, 'n_words': 599, 'path': 'tasks/task_09_CSG_surface_tally_dose/2_surface_dose_from_gamma_source_study.ipynb', 'repo_name': 'pshriwise/openmc_workshop', 'size': '7436'}