added test and improved numerics

This commit is contained in:
GuySten 2026-03-10 23:34:20 +02:00
parent e43e6e7bca
commit 86a7b8f0ce
2 changed files with 61 additions and 8 deletions

View file

@ -669,9 +669,9 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
break;
case SCORE_MIGRATION:
score =
1.0 / 6.0 * p.wgt() *
((p.r() - p.r_born()).norm_sq() - (p.r_last() - p.r_born()).norm_sq());
score = 1.0 / 6.0 * p.wgt() *
(p.r() - p.r_last())
.dot((p.r() - p.r_born()) + (p.r_last() - p.r_born()));
break;
case SCORE_NU_FISSION:
@ -1772,9 +1772,9 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
break;
case SCORE_MIGRATION:
score =
1.0 / 6.0 * p.wgt() *
((p.r() - p.r_born()).norm_sq() - (p.r_last() - p.r_born()).norm_sq());
score = 1.0 / 6.0 * p.wgt() *
(p.r() - p.r_last())
.dot((p.r() - p.r_born()) + (p.r_last() - p.r_born()));
break;
case SCORE_NU_SCATTER:
@ -2646,7 +2646,7 @@ void score_surface_tally(Particle& p, const vector<int>& tallies)
// for a further scoring function.
double score = current * filter_weight;
for (auto score_index = 0; score_index < tally.scores_.size();
++score_index) {
++score_index) {
#pragma omp atomic
tally.results_(filter_index, score_index, TallyResult::VALUE) += score;
}
@ -2721,7 +2721,7 @@ void score_pulse_height_tally(Particle& p, const vector<int>& tallies)
// Loop over scores.
for (auto score_index = 0; score_index < tally.scores_.size();
++score_index) {
++score_index) {
#pragma omp atomic
tally.results_(filter_index, score_index, TallyResult::VALUE) +=
filter_weight;

View file

@ -0,0 +1,53 @@
import openmc
import pytest
import numpy as np
from tests import cdtemp
def sphere_model(radius, boundary_type):
model = openmc.Model()
# Material
material = openmc.Material(name="Hydrogen")
material.add_nuclide("H1", 1.0)
material.set_density('g/cm3', 1.0)
material.add_s_alpha_beta('c_H_in_H2O')
# Geometry
sphere = openmc.Sphere(r=radius, boundary_type=boundary_type)
cell = openmc.Cell(region=-sphere, fill=material)
model.geometry = openmc.Geometry([cell])
# Settings
model.settings.particles = 1000
model.settings.batches = 5
model.settings.run_mode = 'fixed source'
model.settings.source = openmc.IndependentSource()
# Tally
tally = openmc.Tally()
tally.scores = ["migration-area"]
model.tallies = [tally]
return model, tally
def test_reflection_is_equivalent_to_large_model(run_in_tmpdir):
openmc.reset_auto_ids()
refl_model, tally1 = sphere_model(2.5, "reflective")
large_model, tally2 = sphere_model(100, "vacuum")
with cdtemp():
refl_model.run(apply_tally_results=True)
mean1, std1 = (tally1.mean.squeeze(),
tally1.std_dev.squeeze())
with cdtemp():
large_model.run(apply_tally_results=True)
mean2, std2 = (tally2.mean.squeeze(),
tally2.std_dev.squeeze())
std = (std1**2+std2**2)**0.5
diff = np.abs(mean2-mean1)
assert diff <= 3*std