Big difference between MCNP and OpenMC in neutron-photon coupled transport

[YuanHu-PKU-KIT]

Model: 30 cm sphere Material: natural Fe with density 1g/cm3 Source: 55 MeV isotropic neutron at center Tally: photon cell flux Database: for neutron, FENDL-3 was used for both OpenMC and MCNP; for photon, ENDF/B-VII.1 for OpenMC (download by openmc-get-photon-data script) and ENDF/B-VI for MCNP

The results of OpenMC are generally higher than MCNP, with a maximun ratio up to 1.6. But for photon source only, there is a better aggrement between OpenMC and MCNP (within 7%). Does anyone know where the problem might be?

[paulromano]

There are some differences with respect to thick-target bremsstrahlung in the two codes, so it's possible that could be part of the issue. I would try turning off TTB in both codes to see how that affects the results. In OpenMC, you can do this with:

settings = openmc.Settings()
settings.electron_treatment = 'led'
...

[makeclean]

Also at 55 MeV those secondary photons are possible to induce photonuclear, I would also recommend using a dbcn card in mcnp to turn off photonuclear reactions, just in case.

[YuanHu-PKU-KIT]

thick-target bremsstrahlung

I try to turn off TTB in both codes, but the ratio change slightly above 10 MeV, there are still some big difference.

[makeclean]

could you attach a plot? and your inputs - perhaps this discussion should move to the google group?

[YuanHu-PKU-KIT]

Also at 55 MeV those secondary photons are possible to induce photonuclear, I would also recommend using a dbcn card in mcnp to turn off photonuclear reactions, just in case.

The default of PHYS:P card in MCNP is not to include photonuclear collisions, and I also checked the output file, which shows that there are no photons created from photonuclear reactions.

[YuanHu-PKU-KIT]

could you attach a plot? and your inputs - perhaps this discussion should move to the google group?

Here are the OpenMC/MCNP ratio in TTB on and off.

Here are the inputs of OpenMC and MCNP. OpenMC-input.txt MCNP-input.txt

[YuanHu-PKU-KIT]

I used a thin cylinder model (R=0.1μm, 2m longth) to check the results with only single collision, and I found that there are some big differences around 20 MeV.

When I set neutron energy at 19.9999 ~ 19.999999 MeV range, the OpenMC/MCNP ratio are showed below: The OpenMC results are slightly higher than MCNP between 1 ~ 10MeV.

And when I set neutron energy at 19.999999 ~ 20 MeV range, the OpenMC/MCNP ratio are showed below: There are some big difference between OpenMC and MCNP's results. I think this is the reason of the difference I found before. So is this a bug that needed to be fix? And for other elements, does this difference also exist?

P.S: The FENDL-3 database has been used for both MCNP and OpenMC (converted by convert_fendl.py).

[makeclean]

As a general point it would be useful to see which of the variations are statistially important, so in future graphs please include the statistica error too. I wonder if this is a sampling problem, so (correct me if im wrong) I am assuming that you are effectively trying to model a mono-energetic source but your attached source files show linear distributions. This distribution needs to be input per ev which has bitten me before, MCNP does this per bin and does not include the per ev part, which maybe is the reason here? I would suggest trying a neutron source, for mcnp like

sdef x=0 y=0 z=0 erg=20.0

and openmc like

<source>
   <space type="point" x="0" y="0" z="0"/>
  <energy type="discrete" parameters="20e6 1"/>
</source>

This should remove the ambiguity in the source, and perhaps point to a sampling problem.

[YuanHu-PKU-KIT]

As a general point it would be useful to see which of the variations are statistially important, so in future graphs please include the statistica error too. I wonder if this is a sampling problem, so (correct me if im wrong) I am assuming that you are effectively trying to model a mono-energetic source but your attached source files show linear distributions. This distribution needs to be input per ev which has bitten me before, MCNP does this per bin and does not include the per ev part, which maybe is the reason here? I would suggest trying a neutron source, for mcnp like

Yes, I'm trying to model a mono-energetic source. I have changed the input files like you said and run those files again. With the cylinder model, the results are showed bellow:

where red line is the 1 sigma line of MCNP, the blue points are the OpenMC/MCNP ratios. I simulated 1e9 histories, so the statistica error is quite small. And there is still some big difference.

I also repainted those two graphs above:

From 19.99 ~ 19.9999 MeV, the result seems fine at 1 ~ 10 MeV range. But at 19.9999 ~ 20 MeV, there are some difference between OpenMC and MCNP. I also simulated 1e9 histories, so the statistica error is quite small.

[YuanHu-PKU-KIT]

I checked the FENDL3 cross section file for Fe56, the cross section of (n,g) reaction (MT=102) is 0 above 20MeV. So maybe MCNP and OpenMC have different way to deal with this case, which cause this difference.

[makeclean]

I think mcnp uses the last value and leaves it flat, eg if the last value was 20 barns at 20 MeV it will use that for everything up fro

---

From: YuanHu-PKU-KIT notifications@github.com Sent: Thursday, February 20, 2020 3:42:58 PM To: openmc-dev/openmc openmc@noreply.github.com Cc: Davis, Andrew andrew.davis@ukaea.uk; Comment comment@noreply.github.com Subject: Re: [openmc-dev/openmc] Big difference between MCNP and OpenMC in neutron-photon coupled transport (#1472)

I checked the FENDL3 cross section file for Fe56, the cross section of (n,g) reaction (MT=102) is 0 above 20MeV. So maybe MCNP and OpenMC have different way to deal with this case, which cause this difference.

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/openmc-dev/openmc/issues/1472?email_source=notifications&email_token=AASTUSULEXM6Z5UEYRWOXDTRD2QIFA5CNFSM4KQ5Y4K2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEMO2BVY#issuecomment-589144279, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AASTUSQZQZAKDCVIFGZLY53RD2QIFANCNFSM4KQ5Y4KQ.

[YuanHu-PKU-KIT]

I think mcnp uses the last value and leaves it flat, eg if the last value was 20 barns at 20 MeV it will use that for everything up fro … ____ From: YuanHu-PKU-KIT notifications@github.com Sent: Thursday, February 20, 2020 3:42:58 PM To: openmc-dev/openmc openmc@noreply.github.com Cc: Davis, Andrew andrew.davis@ukaea.uk; Comment comment@noreply.github.com Subject: Re: [openmc-dev/openmc] Big difference between MCNP and OpenMC in neutron-photon coupled transport (#1472) I checked the FENDL3 cross section file for Fe56, the cross section of (n,g) reaction (MT=102) is 0 above 20MeV. So maybe MCNP and OpenMC have different way to deal with this case, which cause this difference. — You are receiving this because you commented. Reply to this email directly, view it on GitHub<#1472?email_source=notifications&email_token=AASTUSULEXM6Z5UEYRWOXDTRD2QIFA5CNFSM4KQ5Y4K2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEMO2BVY#issuecomment-589144279>, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AASTUSQZQZAKDCVIFGZLY53RD2QIFANCNFSM4KQ5Y4KQ.

For neutron below 19.9MeV, the n-p coupled transport simulation for OpenMC and MCNP have good agreement, so I may use other database which contain (n,γ) cross sectioin over 20 MeV and to see whether the result will be better.

[paulromano]

Looking at the FENDL 3.1 data, they have cross sections up to 150 MeV so there should be no ambiguity about what happens above 20 MeV. The 0 cross section above 20 MeV should be used by both codes.

@YuanHu-PKU-KIT Rather than using natural Fe, can you try this model with each of the individual isotopes of Fe to see if one in particular is problematic? That could help us narrow down a potential problem.

[YuanHu-PKU-KIT]

When OpenMC samples secondary photons (sample_secondary_photons function in physics.cpp), the number of photons produced y will be decided by photon weight _yt, and then OpenMC will produced y photons which weights are all equal to 1，

//Sample the number of photons produced double y_t = p->wgt_ * p->neutron_xs_[i_nuclide].photon_prod / p->neutron_xs_[i_nuclide].total; int y = static_cast<int>(y_t); if (prn(p->current_seed()) <= y_t - y) ++y;

I modified this code so that each neutron will produce at least 1 photon with weight equal to _yt:

//Sample the number of photons produced double y_t = p->wgt_ * p->neutron_xs_[i_nuclide].photon_prod / p->neutron_xs_[i_nuclide].total; int y = static_cast<int>(y_t); // make sure the secondary photons weight is less than 1, and at least 1 photon will be created y=y+1; double weight=y_t/y; Then set weight to the produced photon. After this modification, the result seems better than before: where the red line is the 3 sigma line of MCNP result, the point is the OpenMC/MCNP ratio. This result is from the same sphere model as the beginning of this issue:

Model: 30 cm sphere Material: natural Fe with density 1g/cm3 Source: 55 MeV isotropic neutron at center Tally: photon cell flux Database: for neutron, FENDL-3 was used for both OpenMC and MCNP; for photon, ENDF/B-VII.1 for OpenMC (download by openmc-get-photon-data script) and ENDF/B-VI for MCNP

But for 20 MeV thin cylinder model, the result is almost the same as before.

[YuanHu-PKU-KIT]

I still could not figure out why this modification works. Those two way should have the same weight in average, but why there are so big difference between two results?

[YuanHu-PKU-KIT]

I found that OpenMC counts the neutron weight twice when created secondary photons, once in the secondary photon weight and once in the number of photons produced. During the OpenMC simulation, the neutron weight is changed in the inelastic scatter process. For neutron energies below 19.9MeV, the neutron inelastic scatter yield is an integer, which causes exactly the same number of neutrons being created. But this yield may be not integer for neutron energies above 19.9MeV, where the reaction MT=5 shows up. At high energies, there are typically many reaction channels open, and it is difficult to decompose the cross section into simple reactions. In such cases, MT=5 reaction has been defined as the sum of all neutron producing reactions not given explicitly elsewhere in the evaluation. For this reaction, OpenMC will assign the non-integer yield as weight to the outgoing neutron. However, when the outgoing neutron with non-integer weight Wn creates a secondary photon in a neutron collision, the expected number of photons produced Np is Np=Wn*σγ/σT where Wn is the neutron weight, σγ is the photon production cross section for the sampled nuclide and σT is the total cross section for the nuclide. [Np] (the truncated integer value of Np) photons will be created with an additional photon produced with probability Np - [Np] and each photon will be assigned the photon weight Wp equal to the neutron weight Wn. So finally the neutron weight has been counted twice. For a neutron weight 1, there is no impact on the results. But for neutron with non-integer weight, this will cause the observed discrepancy. The easiest way to fix this is just not counting the neutron weight Wn through the number of secondary photons as shown in the above formula when photons are generated.

[paulromano]

Nice find @YuanHu-PKU-KIT! I've just confirmed that you can also see this issue below 20 MeV if you turn on survival biasing (settings.survival_biasing = True). Your proposed fixed helps in that situation too.