Absorbed dose rate discrepancy between OpenMC & TRIPOLI

User Support
1

Hi,

I am trying to compare the absorbed dose rate on a cylinder due to a punctual source with both neutron and photon emissions between 3 Monte-Carlo codes: TRIPOLI-4, RayXpert and OpenMC.

Both TRIPOLI-4 and RayXpert gives the same results (around 10% discrepancy, mostly due to some simplification in the TRIPOLI modelling).

However, my OpemnMC input should be identical to the TRIPOLI one but I have a discrepancy factor of 10 000.

Here are the results:

image
image1282×186 7.6 KB

Here is my python script for OpenMC:

#!/opt/conda/bin/python3.10
"""Contain test"""

import datetime

import math
import matplotlib.pyplot as plt
from matplotlib import colors

import openmc

import sys

pwd = os.getcwd()

# Interpolation type
inter_lin_lin = "linear-linear"
inter_log_lin = "log-linear"
inter_lin_log = "linear-log"
inter_log_log = "log-log"

if __name__ == "__main__":#

    # Temperature in K
    t = 293.6

    # Define materials
    # ----------------
    mat = []
    mat_moteur = openmc.Material(material_id=1, name="material_tally", temperature=t)
    mat_moteur.add_element('C',   0.015 , 'wo') 
    mat_moteur.add_element('Si',  0.5   , 'wo')
    mat_moteur.add_element('P',   0.0225, 'wo')
    mat_moteur.add_element('S',   0.0150, 'wo')
    mat_moteur.add_element('Cr',  8.5   , 'wo')
    mat_moteur.add_element('Mn',  1     , 'wo')
    mat_moteur.add_element('Fe', 32.6975, 'wo')
    mat_moteur.add_element('Ni',  6     , 'wo')
    mat_moteur.add_element('Mo',  1.25  , 'wo')
    mat_moteur.add_element('Cu', 50     , 'wo')
    mat_moteur.set_density("g/cm3", 8.45283)
    mat.append(mat_moteur)

    mat_air = openmc.Material(material_id=2, name="Air", temperature=t)
    mat_air.add_element('N', 80, 'ao')
    mat_air.add_element('O', 20, 'ao')
    mat_air.set_density("g/cm3", 1.29E-3)
    mat.append(mat_air)

    # Materials collection
    materials = openmc.Materials(mat)
    
    init=os.getcwd()
    
    for s in ["gamma", "neutron"]:
        os.chdir(init)
        if not os.path.isdir(s):
            os.mkdir(s)
        os.chdir(init+"/"+s)
        materials.export_to_xml()

        # Define geometry = univers = cellules = materials = surfaces
        # -----------------------------------------------------------

        # Surfaces creation
        # moteur
        surf_moteur_ext = openmc.ZCylinder(r=20, x0=-220)
        surf_moteur_zmin = openmc.ZPlane(z0=-20)
        surf_moteur_zmax = openmc.ZPlane(z0=+20)

        # Air ext
        surf_air_ext = openmc.XCylinder(r=50, boundary_type='vacuum')
        surf_air_xmin = openmc.XPlane(x0=-250, boundary_type='vacuum')
        surf_air_xmax = openmc.XPlane(x0=+250, boundary_type='vacuum')

        # Regions creation
        region_moteur = -surf_moteur_ext & -surf_moteur_zmax & +surf_moteur_zmin
        region_air = -surf_air_ext & -surf_air_xmax & +surf_air_xmin
            
        # Cell creation = regions + materials
        cell_tally = openmc.Cell(fill=mat_moteur, region=region_moteur)
        cell_air = openmc.Cell(fill=mat_air, region=region_air)

        # Universe creation
        root_universe = openmc.Universe(universe_id = 1,
                                   name = "Univers",
                                   cells=[cell_tally, cell_air])
        geometry = openmc.Geometry(root_universe)

        # Export
        geometry.export_to_xml()

        # Define settings
        # ---------------
        settings = openmc.Settings()
        # Run mode
        settings.run_mode = 'fixed source'
        # Run strategy
        settings.particles = 20000

        if s == "gamma":
            settings.batches = 100000
        elif s == "neutron":
            settings.batches = 100000

        settings.create_delayed_neutrons = False
        settings.create_fission_neutrons = False

        settings.temperature = {"method": "interpolation"} 
        settings.output = {'tallies': True}

        # Sources
        position = [220,0,0]
        source = openmc.IndependentSource()
        source.space = openmc.stats.Point(position)
        source.angle = openmc.stats.Isotropic()
        
        if s == "neutron":
            # Source neutron
            source.particle = "neutron"
                
            energy_bound = [    
            1.500E-11,
            3.000E-09,
            7.500E-09,
            1.000E-08,
            2.530E-08,
            3.000E-08,
            4.000E-08,
            5.000E-08,
            7.000E-08,
            1.000E-07,
            1.500E-07,
            2.000E-07,
            2.250E-07,
            2.500E-07,
            2.750E-07,
            3.250E-07,
            3.500E-07,
            3.750E-07,
            4.000E-07,
            6.250E-07,
            1.000E-06,
            1.770E-06,
            3.000E-06,
            4.750E-06,
            6.000E-06,
            8.100E-06,
            1.000E-05,
            3.000E-05,
            1.000E-04,
            5.500E-04,
            3.000E-03,
            1.700E-02,
            2.500E-02,
            1.000E-01,
            4.000E-01,
            9.000E-01,
            1.400E+00,
            1.850E+00,
            2.350E+00,
            2.480E+00,
            3.000E+00,
            4.800E+00,
            6.430E+00,
            8.190E+00,
            2.000E+01]
            
            probability = [
            1.328E+08,
            1.048E+09,
            1.233E+09,
            1.383E+10,
            6.686E+09,
            1.772E+10,
            2.203E+10,
            5.285E+10,
            9.419E+10,
            1.725E+11,
            1.707E+11,
            8.055E+10,
            7.485E+10,
            6.893E+10,
            1.211E+11,
            5.598E+10,
            5.176E+10,
            4.740E+10,
            3.051E+11,
            2.944E+11,
            3.450E+11,
            3.178E+11,
            2.765E+11,
            1.395E+11,
            1.806E+11,
            1.280E+11,
            6.580E+11,
            6.808E+11,
            8.461E+11,
            6.766E+11,
            4.960E+11,
            8.628E+10,
            2.468E+11,
            1.227E+11,
            5.755E+10,
            3.067E+10,
            1.382E+10,
            8.411E+09,
            1.412E+09,
            3.225E+09,
            5.480E+09,
            2.693E+09,
            1.235E+09,
            2.884E+08 ]
        elif s == "gamma":
            # Source gamma
            source.particle = "photon"
            
            energy_bound = [    
            1.000E-02,
            2.000E-02,
            5.000E-02,
            7.000E-02,
            1.000E-01,
            1.500E-01,
            2.000E-01,
            2.500E-01,
            3.000E-01,
            3.500E-01,
            4.000E-01,
            4.200E-01,
            4.400E-01,
            4.600E-01,
            4.800E-01,
            5.000E-01,
            5.200E-01,
            5.400E-01,
            5.600E-01,
            5.800E-01,
            6.000E-01,
            6.200E-01,
            6.400E-01,
            6.600E-01,
            6.800E-01,
            7.000E-01,
            7.200E-01,
            7.400E-01,
            7.600E-01,
            7.800E-01,
            8.000E-01,
            8.200E-01,
            8.400E-01,
            8.600E-01,
            8.800E-01,
            8.900E-01,
            9.000E-01,
            9.100E-01,
            9.200E-01,
            9.300E-01,
            9.400E-01,
            9.500E-01,
            9.600E-01,
            9.700E-01,
            9.800E-01,  
            9.900E-01,
            1.000E+00,
            1.100E+00,
            1.200E+00,
            1.300E+00,
            1.400E+00,
            1.500E+00,
            1.600E+00,
            1.700E+00,
            1.800E+00,
            1.900E+00,
            2.000E+00,
            2.100E+00,
            2.200E+00,
            2.300E+00,
            2.400E+00,
            2.500E+00,
            2.600E+00,
            2.700E+00,
            2.800E+00,
            2.900E+00,
            3.000E+00,
            3.200E+00,
            3.400E+00,
            3.600E+00,
            3.800E+00,
            4.000E+00,
            4.500E+00,
            5.000E+00,
            5.500E+00,
            6.000E+00,
            6.500E+00,
            7.000E+00,
            8.000E+00,
            1.000E+01,
            1.200E+01,
            1.400E+01,
            1.600E+01,
            1.800E+01,
            2.000E+01 ]                      
            
            probability = [
            1.389E+08,
            6.037E+09,
            4.089E+10,
            1.130E+11,
            1.723E+11,
            1.265E+11,
            9.420E+10,
            7.663E+10,
            6.537E+10,
            5.945E+10,
            2.207E+10,
            2.058E+10,
            2.033E+10,
            1.059E+11,
            9.525E+09,
            2.883E+10,
            7.594E+09,
            7.305E+09,
            6.726E+09,
            6.632E+09,
            6.444E+09,
            6.311E+09,
            6.516E+09,
            6.202E+09,
            5.871E+09,
            5.273E+09,
            5.314E+09,
            6.601E+09,
            5.670E+09,
            4.799E+09,
            4.422E+09,
            4.939E+09,
            4.921E+09,
            4.258E+09,
            2.230E+09,
            2.462E+09,
            2.206E+09,
            2.096E+09,
            1.922E+09,
            2.242E+09,
            2.083E+09,
            2.076E+09,
            2.340E+09,
            2.209E+09,
            2.223E+09,
            1.957E+09,
            1.875E+10,
            1.744E+10,
            2.446E+10,
            1.471E+10,
            1.358E+10,
            1.267E+10,
            1.184E+10,
            1.611E+10,
            1.070E+10,
            9.505E+09,
            2.304E+10,
            8.271E+09,
            1.050E+10,
            7.997E+09,
            1.034E+10,
            7.515E+09,
            7.283E+09,
            1.019E+10,
            6.688E+09,
            6.309E+09,
            1.664E+10,
            1.180E+10,
            8.665E+10,
            1.320E+10,
            1.077E+10,
            1.458E+10,
            1.051E+11,
            1.032E+10,
            3.881E+09,
            2.197E+10,
            5.117E+09,
            1.565E+10,
            6.740E+09,
            6.825E+08,
            0.000E+00,
            0.000E+00,
            0.000E+00,
            0.000E+00 ]
            
        nb_particle = sum(probability)
        
        energy_bound = [1E6*e for e in energy_bound] # transformation from MeV to eV
        delta_energy = []
        for i in range(len(energy_bound)-1):
            de = energy_bound[i+1]-energy_bound[i]
            delta_energy.append(de)
        
        probability_normalized = [p/e for p, e in zip(probability, delta_energy)]
        probability_normalized.append(0) # adding final 0
        source.energy = openmc.stats.Tabular(x=energy_bound, p=probability_normalized, interpolation="histogram")
        
        settings.source = [source]
        settings.photon_transport = True  # for secondary photon

        # Define tallies
        # --------------
        # Tallies definition
        tallies = openmc.Tallies()
        for p in ["neutron", "electron", "photon", "positron"]:
            tally = openmc.Tally()
            tally.filters = [openmc.CellFilter(bins=cell_tally), openmc.ParticleFilter(p)]
            tally.scores = ["heating"]
            tallies.append(tally)
            
        tally = openmc.Tally()
        tally.filters = [openmc.CellFilter(bins=cell_tally)]
        tally.scores = ["heating"]
        tallies.append(tally)   

        # Export
        tallies.export_to_xml()

        # Export
        settings.export_to_xml()

        # Define model = geometry + materials + settings + tallies
        #--------------
       model = openmc.model.Model(geometry=geometry,
                                   materials=materials,
                                   settings=settings,
                                   tallies=tallies)
       
        openmc.run()

        if 1:
            print("\n########## TALLIES RESULTS ##########")
            with open("resu.txt", "w") as f_out:
                f_out.write(f"Number of particles: {'{:.3E}'.format(nb_particle)}\n")
                f_mult = nb_particle*72800*1.60218e-19/(math.pi*20**2*40*mat_moteur.density/1000)
                f_out.write(f"Multiplication factor is: {'{:.6E}'.format(f_mult)}\n") # from eV/s/source to Gy/s
                               
                with open("tallies.out", "r") as f_in:
                    f_out.write("Résultats : \n")
                    for line in f_in:
                        if "Particle" in line:
                            name = "Particle"
                            f_out.write(f"{'{:10}'.format(line.split()[1].capitalize())} : ")
                        elif "heating" in line:
                            energy_depose = 3600*float(line.split()[1]) # units : from eV/s/source to eV/h/source
                            std_dev = float(line.split()[3])
                            std_dev_pc = 3600*std_dev/energy_depose*100 if energy_depose > 0 else 100.0                                        
                            f_out.write(f"{'{:.3E}'.format(energy_depose)} eV/h/source -> {'{:.3E}'.format(energy_depose*f_mult)} Gy/h +/- {'{:.2f}'.format(std_dev_pc)} %\n") # results in Gy/h

I cannot see where I could be wrong. Any idea is welcome.

2

Bonjour Naja,

you seem to be off by a factor of 10x.

It looks like you are calculating Kerma photon heating (no applied dose function, just absorbed dose).

You should not compare electron and positron heating of OpenMC to the others.
Rather, sum the results of the three (electron heating + photon heating + positron heating) in OpenMC and compare that to Photon Kerma heating in Tripoli.
Electron transport should be off in Tripoli, and ideally both should have similar TTB options. Though that likely doesn’t affect it too much.

I can’t fully follow your normalization.
You appear to multiply the results by the number of particles run? That is atypical.

3

Hi yrrepy,

Thank you for your reply.

Regarding my normalization,

  • nb_particle is the actual number of particles in the neutron or photon spectrum since it’s my understanding that source spectrum is normalized to source.strength (of 1 by default) (I could have put this factor directly in source.strength instead) ;
  • 72800 is a multiplication factor due to the simplification of using a punctual source instead of a surface source (it’s the surface area (in cm2) of the source before simplification ) : it’s the same assumption used in TRIPOLI-4 and the same factor is applied ;
  • the rest is simply transforming the tally results from eV to Gy.

As you mentionned, I am only interested in the total absorbed dose rate (but coming mostly from electron energy deposition).

The version of TRIPOLI-4 used is quite old (version 4.8) and no TTB model exists in this version. Electron and positron should be transported and Bremsstrahlung photon created (according to my knowledge).
However calculations with RayXpert used a TTB model.

Could you explain why TRIPOLI results for electron would be off ,since another french MC code give the same results for electron ?

4

  • So, definitely scrap nb_particle in your normalization
    The result in OpenMC is already normalized to the number of MC particles. You are in fact de-normalizing it.
    I am presuming that you do not have this nb_particle normalization in Tripoli (which it also should not have nor need).

  • The 72800 should then be your source strength (neutrons/second or photons/second).

  • Disable OpenMC TTB
    settings.electron_treatment = 'led'

  • Do not transport electrons and positrons in Tripoli, disable this. OpenMC is not transporting these, so this will not be a fair comparison. The electron, photon and positron bins in the OpenMC tallies is just an accounting of the interaction that lead to the score (tally result). OpenMC does not transport them. The result of their sum should be compared to a Photon heating tally from Tripoli (without electron/positron transport).

  • I mis-added when comparing the tallies, indeed you are off by 10,000x. It is quite likely the nb_particle of the run.