OpenMC_photon_flux_comparison

###############
### Imports ###
###############

import os
import openmc
import numpy as np


#################
### Functions ###
#################

# General function for creating a cuboid region.
def create_cuboid(d, c=[0, 0, 0], r=[0, 0, 0], boundary = "transmission"):
    x_coord, y_coord, z_coord = c # shifting the region
    x_thick, y_thick, z_thick = d # region thicknesses

    # Object rotation
    phi, theta, psi = np.array(r) * np.pi/180
    rotation_transformation = np.round(np.array([
                                [np.cos(theta) * np.cos(psi),
                                -np.cos(phi) * np.sin(psi) + np.sin(phi) * np.sin(theta) * np.cos(psi),
                                 np.sin(phi) * np.sin(psi) + np.cos(phi) * np.sin(theta) * np.cos(psi)],

                                [np.cos(theta) * np.sin(psi),
                                 np.cos(phi) * np.cos(psi) + np.sin(phi) * np.sin(theta) * np.sin(psi),
                                -np.sin(phi) * np.cos(psi) + np.cos(phi) * np.sin(theta) * np.sin(psi)],

                               [-np.sin(theta),
                                 np.sin(phi) * np.cos(theta),
                                 np.cos(phi) * np.cos(theta)]]), 8)

    x_plane_vector = np.array([1, 0, 0])
    x_plane_vector = np.matmul(rotation_transformation, x_plane_vector)

    y_plane_vector = np.array([0, 1, 0])
    y_plane_vector = np.matmul(rotation_transformation, y_plane_vector)

    z_plane_vector = np.array([0, 0, 1])
    z_plane_vector = np.matmul(rotation_transformation, z_plane_vector)

    # Construction of Surface!
    xplus_rect_surface  = openmc.Plane(a=x_plane_vector[0], b=x_plane_vector[1], c=x_plane_vector[2], d=x_thick/2+x_coord, boundary_type=boundary)
    xminus_rect_surface = openmc.Plane(a=x_plane_vector[0], b=x_plane_vector[1], c=x_plane_vector[2], d=-x_thick/2+x_coord, boundary_type=boundary)

    yplus_rect_surface  = openmc.Plane(a=y_plane_vector[0], b=y_plane_vector[1], c=y_plane_vector[2], d=y_thick/2+y_coord, boundary_type=boundary)
    yminus_rect_surface = openmc.Plane(a=y_plane_vector[0], b=y_plane_vector[1], c=y_plane_vector[2], d=-y_thick/2+y_coord, boundary_type=boundary)

    zplus_rect_surface  = openmc.Plane(a=z_plane_vector[0], b=z_plane_vector[1], c=z_plane_vector[2], d=z_thick/2+z_coord, boundary_type=boundary)
    zminus_rect_surface = openmc.Plane(a=z_plane_vector[0], b=z_plane_vector[1], c=z_plane_vector[2], d=-z_thick/2+z_coord, boundary_type=boundary)

    # Sum region
    region = (+xminus_rect_surface & -xplus_rect_surface &
              +yminus_rect_surface & -yplus_rect_surface &
              +zminus_rect_surface & -zplus_rect_surface)

    return region

# User-defined material library
def Material_Library(element_name):
    if element_name == "air":
        mat = openmc.Material(name='air')
        mat.add_element('C', 0.0124, percent_type='wo')
        mat.add_element('N', 75.5267, percent_type='wo')
        mat.add_element('O', 23.1781, percent_type='wo')
        mat.add_element('Ar', 1.2827, percent_type='wo')
        mat.set_density('g/cm3', 1.20479e-3)

    if element_name == "concrete_ordinary":
        mat = openmc.Material(name='concrete_ordinary')
        mat.add_element('H', 1, percent_type='wo')
        mat.add_element('C', 0.1, percent_type='wo')
        mat.add_element('O', 52.9107, percent_type='wo')
        mat.add_element('Na', 1.6, percent_type='wo')
        mat.add_element('Mg', 0.2, percent_type='wo')
        mat.add_element('Al', 3.3872, percent_type='wo')
        mat.add_element('Si', 33.7021, percent_type='wo')
        mat.add_element('K', 1.3, percent_type='wo')
        mat.add_element('Ca', 4.4, percent_type='wo')
        mat.add_element('Fe', 1.4, percent_type='wo')
        mat.set_density('g/cm3', 2.3)

    return mat


#################
### Materials ###
#################

air_mat = Material_Library("air")
concrete_mat = Material_Library("concrete_ordinary")
materials = openmc.Materials([air_mat, concrete_mat])

materials.export_to_xml()


################
### Settings ###
################

# Run Settings
settings = openmc.Settings()
settings.batches = 10
settings.inactive = 0
settings.particles = 50000000
settings.run_mode = "fixed source"

# Transport Settings
settings.max_tracks = 1000
settings.electron_treatment = 'ttb' # led ttb

# initialises a new source object.
my_source = openmc.Source()

# the distribution of source z values is just a single value.
x_values = openmc.stats.Uniform(a=-50, b=50)
y_values = openmc.stats.Uniform(a=-50, b=50)
z_values = openmc.stats.Uniform(a=0, b=0)

my_source.space = openmc.stats.CartesianIndependent(x=x_values, y=y_values, z=z_values)

# sets the direction normal to concrete surface.
my_source.angle = openmc.stats.Monodirectional([0, 0, 1])

# sets the energy distribution.
my_source.energy = openmc.stats.Discrete([0.6616e6, 1.1732e6, 1.3325e6], [0.8605, 0.0697, 0.0698])
my_source.particle = 'photon'

# Select Source for OpenMC run.
settings.source = my_source

settings.cutoff = {'energy_electron': 0.521e6,
                   'energy_photon'  : 1e4}

# Export to xml.
settings.export_to_xml()


################
### Geometry ###
################

left_layer_value = 0       # For shifting layers to the correct positions.
concrete_thickness = 10    # thickness of concrete layers in cm.

# Initial air layer of 150 cm.
left_layer_value += 75
ini_air_region = create_cuboid([300, 300, 150], c=[0, 0, left_layer_value])
ini_air_cell = openmc.Cell(region=ini_air_region)
ini_air_cell.fill = air_mat
left_layer_value += 75

# Concrete layers.
left_layer_value += concrete_thickness/2
concrete_region_1 = create_cuboid(d=[300, 300, concrete_thickness], c=[0, 0, left_layer_value])
concrete_cell_1 = openmc.Cell(region=concrete_region_1)
concrete_cell_1.fill = concrete_mat
left_layer_value += concrete_thickness/2

left_layer_value += concrete_thickness/2
concrete_region_2 = create_cuboid(d=[300, 300, concrete_thickness], c=[0, 0, left_layer_value])
concrete_cell_2 = openmc.Cell(region=concrete_region_2)
concrete_cell_2.fill = concrete_mat
left_layer_value += concrete_thickness/2

left_layer_value += concrete_thickness/2
concrete_region_3 = create_cuboid(d=[300, 300, concrete_thickness], c=[0, 0, left_layer_value])
concrete_cell_3 = openmc.Cell(region=concrete_region_3)
concrete_cell_3.fill = concrete_mat
left_layer_value += concrete_thickness/2

left_layer_value += concrete_thickness/2
concrete_region_4 = create_cuboid(d=[300, 300, concrete_thickness], c=[0, 0, left_layer_value])
concrete_cell_4 = openmc.Cell(region=concrete_region_4)
concrete_cell_4.fill = concrete_mat
left_layer_value += concrete_thickness/2

left_layer_value += concrete_thickness/2
concrete_region_5 = create_cuboid(d=[300, 300, concrete_thickness], c=[0, 0, left_layer_value])
concrete_cell_5 = openmc.Cell(region=concrete_region_5)
concrete_cell_5.fill = concrete_mat
left_layer_value += concrete_thickness/2

# Final air layers of 10 cm and 140 cm.
left_layer_value += 5
fin_air_region_1 = create_cuboid([300, 300, 10], c=[0, 0, left_layer_value])
fin_air_cell_1 = openmc.Cell(region=fin_air_region_1)
fin_air_cell_1.fill = air_mat
left_layer_value += 5

left_layer_value += 70
fin_air_region_2 = create_cuboid([300, 300, 140], c=[0, 0, left_layer_value])
fin_air_cell_2 = openmc.Cell(region=fin_air_region_2)
fin_air_cell_2.fill = air_mat
left_layer_value += 70

# Vacuum
outer_region = create_cuboid([300.001, 300.001, left_layer_value+0.001], c=[0, 0, left_layer_value/2], boundary="vacuum")
outer_cell = openmc.Cell(region=outer_region)

# makes a universe to contain all the cells
universe = openmc.Universe(cells=[ini_air_cell,
                                  concrete_cell_1, concrete_cell_2, concrete_cell_3, concrete_cell_4, concrete_cell_5,
                                  fin_air_cell_1, fin_air_cell_2,
                                  outer_cell])
geometry = openmc.Geometry(universe)

geometry.export_to_xml()


###############
### Tallies ###
###############

tallies = openmc.Tallies()

particle_filter = openmc.ParticleFilter("photon")
cell_filter = openmc.CellFilter([concrete_cell_1, concrete_cell_2, concrete_cell_3, concrete_cell_4, concrete_cell_5,fin_air_cell_1])
energy_array = np.array([0.01, 0.0249, 0.0398, 0.0547, 0.0696, 0.0845, 0.0994, 0.1143, 0.1292, 0.1441, 0.159, 0.1739, 0.1888,
                         0.2037, 0.2186, 0.2335, 0.2484, 0.2633, 0.2782, 0.2931, 0.308, 0.3229, 0.3378, 0.3527, 0.3676, 0.3825,
                         0.3974, 0.4123, 0.4272, 0.4421, 0.457, 0.4719, 0.4868, 0.5017, 0.5166, 0.5315, 0.5464, 0.5613, 0.5762,
                         0.5911, 0.606, 0.6209, 0.6358, 0.6507, 0.6656, 0.6805, 0.6954, 0.7103, 0.7252, 0.7401, 0.755, 0.7699,
                         0.7848, 0.7997, 0.8146, 0.8295, 0.8444, 0.8593, 0.8742, 0.8891, 0.904, 0.9189, 0.9338, 0.9487, 0.9636,
                         0.9785, 0.9934, 1.0083, 1.0232, 1.0381, 1.053, 1.0679, 1.0828, 1.0977, 1.1126, 1.1275, 1.1424, 1.1573,
                         1.1722, 1.1871, 1.202, 1.2169, 1.2318, 1.2467, 1.2616, 1.2765, 1.2914, 1.3063, 1.3212, 1.3361, 1.351,
                         1.3659, 1.3808, 1.3957, 1.4106, 1.4255, 1.4404, 1.4553, 1.4702, 1.4851, 1.5])
energy_array *= 1e6
energy_filter = openmc.EnergyFilter(energy_array)

flux_tally = openmc.Tally(name='Photon Flux')
flux_tally.filters = [cell_filter, energy_filter, particle_filter]

flux_tally.scores = ['flux']
tallies.append(flux_tally)


###################
### Calculation ###
###################

model = openmc.model.Model(geometry, materials, settings, tallies)
sp_filename = model.run(tracks=True)


###############
### Results ###
###############

# calculate mean flux per volume
tally = openmc.StatePoint('statepoint.10.h5').get_tally(name="Photon Flux")
photon_flux = tally.mean
photon_flux = photon_flux / (300*300*10)