Dear OpenMC helpers,
I’m new to openmc (0.14.0) and i couldn’t found any similar problem …
I got negative values for neutron flux scores … i can’t find out what is wrong …
the code 
import sys
import json
import time
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import openmc
import openmc.lib
import neutronics_material_maker as nmm
Matériaux
ss316 = openmc.Material(name=‘SS316’)
ss316.set_density(‘g/cm3’, 8.0)
ss316.add_nuclide(‘Fe56’, 0.65)
ss316.add_nuclide(‘Cr52’, 0.17)
ss316.add_nuclide(‘Ni58’, 0.12)
ss316.add_nuclide(‘Mo98’, 0.025)
ss316.add_nuclide(‘Mn55’, 0.02)
ss316.set_density(‘g/cm3’, 7.98)
tungsten = openmc.Material(name=‘Tungsten’)
tungsten.add_element(‘W’, 1.0)
tungsten.set_density(‘g/cm3’, 19.25)
graphite = openmc.Material(name=‘Graphite’)
graphite.add_element(‘C’, 1.0)
graphite.set_density(‘g/cm3’, 1.7)
carbon_steel = openmc.Material(name=‘Carbon Steel’)
carbon_steel.add_element(‘Fe’, 0.98)
carbon_steel.add_element(‘C’, 0.02)
carbon_steel.set_density(‘g/cm3’, 7.85)
tungsten = openmc.Material(name=‘Tungsten’)
tungsten.add_element(‘W’, 1.0)
tungsten.set_density(‘g/cm3’, 19.25)
tritium_source = openmc.Material(name=‘TritiumSource’)
tritium_source.add_nuclide(‘H3’, 1.0)
air = openmc.Material(name=‘Air’)
air.add_element(‘N’, 0.78) # Nitrogen makes up about 78% of air by volume
air.add_element(‘O’, 0.21) # Oxygen makes up about 21% of air by volume
air.add_element(‘Ar’, 0.01) # Argon makes up about 1% of air by volume
air.set_density(‘g/cm3’, 0.001225)
Define natural lithium
natural_li = openmc.Material(name=‘Natural Lithium’)
natural_li.add_nuclide(‘Li6’, 0.0759)
natural_li.add_nuclide(‘Li7’, 0.9241)
natural_li.set_density(‘g/cm3’, 0.534)
materials = openmc.Materials([ss316, tungsten, graphite, carbon_steel, tritium_source, natural_li, air])
Surfaces
chamber_outer = openmc.ZCylinder(x0=0, y0=0, r=60, boundary_type=“vacuum”) # boundary_type=‘reflective’
chamber_inner = openmc.ZCylinder(x0=0, y0=0, r=55)
graphite_region = openmc.ZCylinder(x0=0, y0=0, r=50)
carbon_steel_vessel = openmc.ZCylinder(x0=0, y0=0, r=45)
tungsten_crucible = openmc.ZCylinder(x0=0, y0=0, r=30)
breeder_region = openmc.ZCylinder(x0=0, y0=0, r=20)
source_cavity = openmc.ZCylinder(x0=0, y0=0, r=10)
front_p = openmc.ZPlane(z0=-100/2, name=‘Front plane’ , boundary_type=‘reflective’)
back_p = openmc.ZPlane(z0= 100/2, name=‘Back plane’ , boundary_type=‘reflective’)
surf_Xplus = openmc.XPlane(x0= 0, name=‘separation plane’, boundary_type=‘reflective’)
Cellules
cell1 = openmc.Cell(region= -chamber_outer & +chamber_inner & +front_p & -back_p & +surf_Xplus, fill=air)
cell2 = openmc.Cell(region= -chamber_inner & +graphite_region & +front_p & -back_p & +surf_Xplus, fill=ss316)
cell3 = openmc.Cell(region= -graphite_region & +carbon_steel_vessel & +front_p & -back_p & +surf_Xplus, fill=graphite)
cell4 = openmc.Cell(region= -carbon_steel_vessel & +tungsten_crucible & +front_p & -back_p & +surf_Xplus, fill=carbon_steel)
cell5 = openmc.Cell(region= -tungsten_crucible & +breeder_region & +front_p & -back_p & +surf_Xplus, fill=tungsten)
cell6 = openmc.Cell(region= -breeder_region & +source_cavity & +front_p & -back_p & +surf_Xplus, fill=natural_li)
cell7 = openmc.Cell(region= -source_cavity & +front_p & -back_p & +surf_Xplus, fill=air)
Géométrie
geometry = openmc.Geometry()
Création de l’univers racine
root = openmc.Universe(universe_id=0, name=‘root universe’)
Ajout des cellules à l’univers racine
root.add_cells([cell1, cell2, cell3, cell4, cell5, cell6, cell7])
Création de la géométrie avec l’univers racine
geometry = openmc.Geometry(root)
source = openmc.IndependentSource()
source.particle=‘neutron’
rad_src = openmc.stats.Uniform(a=0, b=10.)
phi_src = openmc.stats.Uniform(a=0, b=23.14159265359) # b=23.14159265359
z_src = openmc.stats.Uniform(a=-10/2, b=10/2)
origin_src = (0.0, 0.0, 0.0)
source.space = openmc.stats.CylindricalIndependent(r=rad_src, phi=phi_src, z=z_src,origin=origin_src)
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Discrete([141e5], [1])
Paramètres de simulation
settings = openmc.Settings()
settings.run_mode = “fixed source”
settings.source = source
settings.batches = 10
settings.particles = 10000
settings.output = {‘tallies’: False}
Tallies
mesh = openmc.CylindricalMesh.from_domain(
domain=geometry, # the corners of the mesh are being set automatically to surround the geometry
dimension=[60, 1, 1] # 60 voxels in r axis direction
)
mesh_filter = openmc.MeshFilter(mesh)
particle_filter = openmc.ParticleFilter(‘neutron’)
cell_filter = openmc.CellFilter([cell1, cell2, cell3, cell4, cell5, cell6, cell7])
tally = openmc.Tally(tally_id=1)
tally.filters = [mesh_filter, cell_filter ]# particle_filter
tally.filters = [mesh_filter ]# particle_filter
tally.scores = [‘flux’]
tallies = openmc.Tallies([tally])
Modèle OpenMC
model = openmc.Model(geometry, materials, settings, tallies)
Exporter au format XML
model.export_to_xml()
deletes old statepoint and summary files
!rm *.h5
Run OpenMC!
model = openmc.model.Model(geometry, materials, settings, tallies)
sp_filename = model.run()
sp = openmc.StatePoint(‘statepoint.10.h5’)
tally = sp.get_tally(scores=[‘flux’])
print(tally)
tally.sum
print(tally.mean.shape)
mean_values = tally.mean.flatten()
std_dev_values = tally.std_dev.flatten()
data = {
‘mean_values’: mean_values,
‘std_dev_values’: std_dev_values
}
df = pd.DataFrame(data)
df1 = tally.get_pandas_dataframe(nuclides=False)
Display the DataFr
ame in a more readable format
print(df)
print(df1.head(10))
df1[df1[“mean_values”] > 0]
Extract the coordinates of the mesh points
r_edges = mesh.r_grid
r_centers = 0.5 * (r_edges[:-1] + r_edges[1:]) # Compute the center points of each radial bin
print( geometry.find( np.array(r_centers[0],0.,0.) ) )
print( geometry.find( (float(r_centers[0]), 0., 0.))[1].id )
Create a plot of distance vs radius
plt.figure(figsize=(8, 6))
plt.plot(r_centers, mean_values, marker=‘o’)
plt.xlabel(‘Radius (R)’)
plt.ylabel(‘Flux [a.u.]’)
plt.title(‘Flux vs Radius (R)’)
plt.grid(True)
plt.show()
resulting
Tally
ID = 1
Name =
Filters = MeshFilter, CellFilter
Nuclides = total
Scores = [‘flux’]
Estimator = tracklength
Multiply dens. = True
(420, 1, 1)
mean_values std_dev_values
0 0.000000 0.000000
1 0.000000 0.000000
2 0.000000 0.000000
3 -4.798954 1.623547
4 -4.390122 0.749598
… … …
415 0.000000 0.000000
416 0.000020 0.000014
417 0.000000 0.000000
418 0.000000 0.000000
419 0.000000 0.000000
[420 rows x 2 columns]
mesh 2 cell score mean std. dev.
x y z
0 1 1 1 8 flux 0.000000 0.000000
1 1 1 1 9 flux 0.000000 0.000000
2 1 1 1 10 flux 0.000000 0.000000
3 1 1 1 11 flux -4.798954 1.623547
4 1 1 1 12 flux -4.390122 0.749598
5 1 1 1 13 flux 0.000000 0.000000
6 1 1 1 14 flux 0.000000 0.000000
7 2 1 1 8 flux 0.000000 0.000000
8 2 1 1 9 flux 0.000000 0.000000
9 2 1 1 10 flux 0.000000 0.000000
14
=> thus negative values for the flux !!!