Photoelectron direction sampled in the lab frame (fixed +x axis) instead of relative to the incident photon — creates a spurious +x bremsstrahlung beam with TTB

オープン 初心者向け
  • 何か問題がありますか？
#4,222 コメント 1 件 リアクション 0 件 担当者 0 名 GitHub で見る

メンテナーはふだん 2 日以内に返信

まだ誰も着手していません。

評価

難易度
2/5
見積もり時間
1〜3時間
初心者へのやさしさ
78/100
issue の種類
バグ
明瞭さ
明確に書かれている
活発さ
活発
技術スタック
cpp, python
領域
backend

調査の方向性

src/physics.cpp の sample_photon_reaction の光電効果分岐を読み、方向の扱いをコンプトン散乱分岐および対生成分岐と比較してください。Issue には、指定された鉛殻の形状と electron_treatment = "ttb" を使った修正案と Python による再現例が示されています。変更後、6つの検出器球が球殻の基準値と一致することを確認してください。

索引モデルが issue の本文から書いたものです。

説明

Bugs

Bug Description

In sample_photon_reaction (src/physics.cpp, photoelectric branch), the photoelectron emission cosine mu is sampled from the Sauter distribution, which is defined relative to the incident photon direction. The resulting direction is however built directly in the lab frame:

double phi = uniform_distribution(0., 2.0 * PI, p.current_seed());
Direction u;
u.x = mu;
u.y = std::sqrt(1.0 - mu * mu) * std::cos(phi);
u.z = std::sqrt(1.0 - mu * mu) * std::sin(phi);

process_charged_secondary(p, u, E_electron, ParticleType::electron());

so every photoelectron is emitted around the fixed +x axis, whatever the photon direction. The Compton and pair-production branches correctly use rotate_angle(p.u(), ...).

With electron_treatment = "led" the electron direction is never used, so there is no visible effect. With "ttb" (the default), thick_target_bremsstrahlung emits the bremsstrahlung photons along the electron direction u, so the photoelectric bremsstrahlung becomes a beam along +x. At high photon energies the Sauter distribution is strongly forward-peaked (mu → 1), so the beam is narrow and energetic. In high-Z shields (Pb, W) at MeV energies this is clearly visible at detectors located on the +x side of the source, while other directions are slightly under-estimated.

Present in v0.16.0 and in develop.

Proposed fix
         double phi = uniform_distribution(0., 2.0 * PI, p.current_seed());
-        Direction u;
-        u.x = mu;
-        u.y = std::sqrt(1.0 - mu * mu) * std::cos(phi);
-        u.z = std::sqrt(1.0 - mu * mu) * std::sin(phi);
+        Direction u = rotate_angle(p.u(), mu, &phi, p.current_seed());

Happy to open a PR if that helps.

Steps to Reproduce

Spherically symmetric problem: isotropic 2614 keV point source at the centre of a lead shell (r = 1 to 6 cm), vacuum elsewhere. Six identical void spheres (r = 3 cm) at 90 cm on ±x, ±y, ±z, plus a thin spherical shell at 65 cm (4π average, scaled by 1/r²). By symmetry all six spheres must agree with each other and with the shell.

openmc_repro_photoelectron.py 
"""Minimal reproducer: photoelectron direction sampled in the lab frame (OpenMC issue draft).

Isotropic 2614 keV point source at the centre of a lead sphere (r = 1 to 6 cm), vacuum outside.
The problem is spherically symmetric, so six identical void spheres (r = 3 cm) at 65 cm on
±x, ±y and ±z must all see the same photon flux, and the same as a thin spherical shell at
65 cm. With electron_treatment = "ttb", the +x sphere reads high.

    python openmc_repro_photoelectron.py [ttb|led] [batches]
"""
import math
import sys
import tempfile

import numpy as np
import openmc

treatment = sys.argv[1] if len(sys.argv) > 1 else "ttb"
batches = int(sys.argv[2]) if len(sys.argv) > 2 else 400

pb = openmc.Material()
pb.add_element("Pb", 1.0)
pb.set_density("g/cm3", 11.35)

inner, outer = openmc.Sphere(r=1.0), openmc.Sphere(r=6.0)
s0, s1 = openmc.Sphere(r=64.5), openmc.Sphere(r=65.5)
world = openmc.Sphere(r=120.0, boundary_type="vacuum")
dirs = [(65, 0, 0), (-65, 0, 0), (0, 65, 0), (0, -65, 0), (0, 0, 65), (0, 0, -65)]
# the balls sit outside the thin shell (centres at 65 cm are inside it), so put them at 90 cm
dirs = [tuple(90 * c / 65 for c in d) for d in dirs]
balls = [openmc.Sphere(x0=x, y0=y, z0=z, r=3.0) for x, y, z in dirs]

ball_cells = [openmc.Cell(region=-b) for b in balls]
gap = +s1 & -world
for b in balls:
    gap &= +b
shell_cell = openmc.Cell(region=+s0 & -s1)
cells = [openmc.Cell(region=-inner), openmc.Cell(fill=pb, region=+inner & -outer),
         openmc.Cell(region=+outer & -s0), shell_cell, *ball_cells, openmc.Cell(region=gap)]

settings = openmc.Settings()
settings.run_mode = "fixed source"
settings.particles = 200_000
settings.batches = batches
settings.photon_transport = True
settings.electron_treatment = treatment
settings.seed = 1
settings.source = openmc.IndependentSource(space=openmc.stats.Point(), angle=openmc.stats.Isotropic(),
                                           energy=openmc.stats.Discrete([2.6145e6], [1.0]), particle="photon")

tallies = openmc.Tallies()
for name, cs in (("balls", ball_cells), ("shell", [shell_cell])):
    t = openmc.Tally(name=name)
    t.filters = [openmc.CellFilter(cs), openmc.ParticleFilter(["photon"])]
    t.scores = ["flux"]
    tallies.append(t)

model = openmc.Model(openmc.Geometry(cells), openmc.Materials([pb]), settings, tallies)
with tempfile.TemporaryDirectory() as d:
    sp_path = model.run(cwd=d, output=False)
    with openmc.StatePoint(sp_path) as sp:
        b, s = sp.get_tally(name="balls"), sp.get_tally(name="shell")
        vb = 4 / 3 * math.pi * 3.0 ** 3
        vs = 4 / 3 * math.pi * (65.5 ** 3 - 64.5 ** 3)
        ref = s.mean.ravel()[0] / vs * (65 / 90) ** 2  # shell flux scaled to 90 cm
        print(f"OpenMC {openmc.__version__}, electron_treatment = {treatment}, {batches} x 200000 photons")
        print(f"thin shell (4π average), scaled to 90 cm: {ref:.4e} ± {100 * s.std_dev.ravel()[0] / s.mean.ravel()[0]:.2f} %")
        for (x, y, z), m, e in zip(dirs, b.mean.ravel(), b.std_dev.ravel()):
            lab = "+x" if x > 0 else "-x" if x < 0 else "+y" if y > 0 else "-y" if y < 0 else "+z" if z > 0 else "-z"
            print(f"sphere {lab}: {m / vb:.4e} ± {100 * e / m:.1f} %   ratio to shell {m / vb / ref:.3f}")
Results — v0.16.0 (conda-forge), ENDF/B-VIII.0 photon data
OpenMC 0.16.0, electron_treatment = ttb, 400 x 200000 photons
thin shell (4π average), scaled to 90 cm: 2.0935e-06 ± 0.02 %
sphere +x: 2.3528e-06 ± 1.4 %   ratio to shell 1.124
sphere -x: 2.0863e-06 ± 1.6 %   ratio to shell 0.997
sphere +y: 2.0684e-06 ± 1.6 %   ratio to shell 0.988
sphere -y: 2.0685e-06 ± 1.5 %   ratio to shell 0.988
sphere +z: 2.1247e-06 ± 1.6 %   ratio to shell 1.015
sphere -z: 2.0889e-06 ± 1.6 %   ratio to shell 0.998
Results — same, with the patch below
OpenMC 0.16.0 + patch, electron_treatment = ttb, 400 x 200000 photons
thin shell (4π average), scaled to 90 cm: 2.1135e-06 ± 0.02 %
sphere +x: 2.0862e-06 ± 1.5 %   ratio to shell 0.987
sphere -x: 2.1266e-06 ± 1.6 %   ratio to shell 1.006
sphere +y: 2.1005e-06 ± 1.6 %   ratio to shell 0.994
sphere -y: 2.1022e-06 ± 1.5 %   ratio to shell 0.995
sphere +z: 2.1549e-06 ± 1.5 %   ratio to shell 1.020
sphere -z: 2.1219e-06 ± 1.6 %   ratio to shell 1.004

We first noticed it in a real shielding problem (Pb-212 / Tl-208 2614 keV behind 5 cm of lead): a point-like detector on +x read ~10 % above the 4π-averaged value, and the excess disappeared with electron_treatment = "led" or with a 662 keV source.

The +x sphere is 12 % (~8σ) above the 4π average before the patch and consistent with it after. Note that the 4π-averaged value itself also changes (+1.0 %): with the bremsstrahlung emitted along the physical photoelectron direction (mostly outward), more of it escapes the shield. So integral results are affected too, not only point-like tallies near the +x axis.

Environment

  • OpenMC 0.16.0 (conda-forge build dagmc_nompi_py314), Linux (WSL2 Ubuntu 24.04)
  • Cross sections: ENDF/B-VIII.0 HDF5 (official OpenMC library)
  • The same code is still present in develop (src/physics.cpp, photoelectric branch of sample_photon_reaction, checked 2026-11-08).

openmc-photoelectron-direction.patch

主要言語
Python
スター
1.1k
フォーク
683
平均マージ
3日 10時間
マージ済み PR（30日）
17

環境構築

はじめの一歩

  1. issue を最後まで読み、次にプロジェクトのコントリビューションガイドを読みます。
  2. 着手することを issue にコメントします — 二人が同じ作業をするのを防げます。
  3. リポジトリをフォークし、ブランチを切って変更します。
  4. issue 番号を参照したプルリクエストを送ります。

openmc-dev/openmc のほかの issue

openmc-dev/openmc の issue をすべて見る

似ている issue

Python の issue をもっと見る

新しい issue をメールで受け取る

初心者向けの GitHub issue を短くまとめたダイジェスト。