crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Fix energy distributions to correctly sample discrete lines

Browse source
This commit is contained in:
amandalund 2018-09-11 08:46:12 -05:00
parent 6db62984d7
commit 3f0af69aef
4 changed files with 102 additions and 49 deletions
Download patch file Download diff file Expand all files Collapse all files

39
src/distribution_energy.cpp
View file

@ -182,11 +182,13 @@ double ContinuousTabular::sample(double E) const
// Interpolation for energy E1 and EK
int n_energy_out = distribution_[i].e_out.size();
double E_i_1 = distribution_[i].e_out[0];
int n_discrete = distribution_[i].n_discrete;
double E_i_1 = distribution_[i].e_out[n_discrete];
double E_i_K = distribution_[i].e_out[n_energy_out - 1];
n_energy_out = distribution_[i+1].e_out.size();
double E_i1_1 = distribution_[i+1].e_out[0];
n_discrete = distribution_[i+1].n_discrete;
double E_i1_1 = distribution_[i+1].e_out[n_discrete];
double E_i1_K = distribution_[i+1].e_out[n_energy_out - 1];
double E_1 = E_i_1 + r*(E_i1_1 - E_i_1);
@ -194,25 +196,38 @@ double ContinuousTabular::sample(double E) const
// Determine outgoing energy bin
n_energy_out = distribution_[l].e_out.size();
n_discrete = distribution_[l].n_discrete;
double r1 = prn();
double c_k = distribution_[l].c[0];
double c_k1;
int k;
for (k = 0; k < n_energy_out - 2; ++k) {
c_k1 = distribution_[l].c[k+1];
if (r1 < c_k1) break;
c_k = c_k1;
int k = 0;
int end = n_energy_out - 2;
// Discrete portion
for (int j = 0; j < n_discrete; ++j) {
k = j;
c_k = distribution_[l].c[k];
if (r1 < c_k) {
end = j;
break;
}
}
// Check to make sure 1 <= k <= NP - 1
k = std::max(0, std::min(k, n_energy_out - 2));
// Continuous portion
double c_k1;
for (int j = n_discrete; j < end; ++j) {
k = j;
c_k1 = distribution_[l].c[k+1];
if (r1 < c_k1) break;
k = j + 1;
c_k = c_k1;
}
double E_l_k = distribution_[l].e_out[k];
double p_l_k = distribution_[l].p[k];
double E_out;
if (distribution_[l].interpolation == Interpolation::histogram) {
// Histogram interpolation
if (p_l_k > 0.0 && k >= distribution_[l].n_discrete) {
if (p_l_k > 0.0 && k >= n_discrete) {
E_out = E_l_k + (r1 - c_k)/p_l_k;
} else {
E_out = E_l_k;
@ -233,7 +248,7 @@ double ContinuousTabular::sample(double E) const
}
// Now interpolate between incident energy bins i and i + 1
if (!histogram_interp && n_energy_out > 1) {
if (!histogram_interp && n_energy_out > 1 && k >= n_discrete) {
if (l == i) {
return E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1);
} else {

18
src/photon_header.F90
View file

@ -163,14 +163,18 @@ contains
call this % coherent_int_form_factor % from_hdf5(dset_id)
call close_dataset(dset_id)
dset_id = open_dataset(rgroup, 'anomalous_real')
call this % coherent_anomalous_real % from_hdf5(dset_id)
call close_dataset(dset_id)
if (object_exists(group_id, 'anomalous_real')) then
dset_id = open_dataset(rgroup, 'anomalous_real')
call this % coherent_anomalous_real % from_hdf5(dset_id)
call close_dataset(dset_id)
end if
dset_id = open_dataset(rgroup, 'anomalous_imag')
call this % coherent_anomalous_imag % from_hdf5(dset_id)
call close_dataset(dset_id)
call close_group(rgroup)
if (object_exists(group_id, 'anomalous_imag')) then
dset_id = open_dataset(rgroup, 'anomalous_imag')
call this % coherent_anomalous_imag % from_hdf5(dset_id)
call close_dataset(dset_id)
call close_group(rgroup)
end if
! Read incoherent scattering
rgroup = open_group(group_id, 'incoherent')

47
src/secondary_correlated.cpp
View file

@ -177,11 +177,13 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu) const
// Interpolation for energy E1 and EK
int n_energy_out = distribution_[i].e_out.size();
double E_i_1 = distribution_[i].e_out[0];
int n_discrete = distribution_[i].n_discrete;
double E_i_1 = distribution_[i].e_out[n_discrete];
double E_i_K = distribution_[i].e_out[n_energy_out - 1];
n_energy_out = distribution_[i+1].e_out.size();
double E_i1_1 = distribution_[i+1].e_out[0];
n_discrete = distribution_[i+1].n_discrete;
double E_i1_1 = distribution_[i+1].e_out[n_discrete];
double E_i1_K = distribution_[i+1].e_out[n_energy_out - 1];
double E_1 = E_i_1 + r*(E_i1_1 - E_i_1);
@ -189,24 +191,37 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu) const
// Determine outgoing energy bin
n_energy_out = distribution_[l].e_out.size();
n_discrete = distribution_[l].n_discrete;
double r1 = prn();
double c_k = distribution_[l].c[0];
double c_k1;
int k;
for (k = 0; k < n_energy_out - 2; ++k) {
c_k1 = distribution_[l].c[k+1];
if (r1 < c_k1) break;
c_k = c_k1;
int k = 0;
int end = n_energy_out - 2;
// Discrete portion
for (int j = 0; j < n_discrete; ++j) {
k = j;
c_k = distribution_[l].c[k];
if (r1 < c_k) {
end = j;
break;
}
}
// Check to make sure 1 <= k <= NP - 1
k = std::max(0, std::min(k, n_energy_out - 2));
// Continuous portion
double c_k1;
for (int j = n_discrete; j < end; ++j) {
k = j;
c_k1 = distribution_[l].c[k+1];
if (r1 < c_k1) break;
k = j + 1;
c_k = c_k1;
}
double E_l_k = distribution_[l].e_out[k];
double p_l_k = distribution_[l].p[k];
if (distribution_[l].interpolation == Interpolation::histogram) {
// Histogram interpolation
if (p_l_k > 0.0 && k >= distribution_[l].n_discrete) {
if (p_l_k > 0.0 && k >= n_discrete) {
E_out = E_l_k + (r1 - c_k)/p_l_k;
} else {
E_out = E_l_k;
@ -228,10 +243,12 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu) const
}
// Now interpolate between incident energy bins i and i + 1
if (l == i) {
E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1);
} else {
E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1);
if (k >= n_discrete){
if (l == i) {
E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1);
} else {
E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1);
}
}
// Find correlated angular distribution for closest outgoing energy bin

47
src/secondary_kalbach.cpp
View file

@ -144,11 +144,13 @@ void KalbachMann::sample(double E_in, double& E_out, double& mu) const
// Interpolation for energy E1 and EK
int n_energy_out = distribution_[i].e_out.size();
double E_i_1 = distribution_[i].e_out[0];
int n_discrete = distribution_[i].n_discrete;
double E_i_1 = distribution_[i].e_out[n_discrete];
double E_i_K = distribution_[i].e_out[n_energy_out - 1];
n_energy_out = distribution_[i+1].e_out.size();
double E_i1_1 = distribution_[i+1].e_out[0];
n_discrete = distribution_[i+1].n_discrete;
double E_i1_1 = distribution_[i+1].e_out[n_discrete];
double E_i1_K = distribution_[i+1].e_out[n_energy_out - 1];
double E_1 = E_i_1 + r*(E_i1_1 - E_i_1);
@ -156,25 +158,38 @@ void KalbachMann::sample(double E_in, double& E_out, double& mu) const
// Determine outgoing energy bin
n_energy_out = distribution_[l].e_out.size();
n_discrete = distribution_[l].n_discrete;
double r1 = prn();
double c_k = distribution_[l].c[0];
double c_k1;
int k;
for (k = 0; k < n_energy_out - 2; ++k) {
c_k1 = distribution_[l].c[k+1];
if (r1 < c_k1) break;
c_k = c_k1;
int k = 0;
int end = n_energy_out - 2;
// Discrete portion
for (int j = 0; j < n_discrete; ++j) {
k = j;
c_k = distribution_[l].c[k];
if (r1 < c_k) {
end = j;
break;
}
}
// Check to make sure 1 <= k <= NP - 1
k = std::max(0, std::min(k, n_energy_out - 2));
// Continuous portion
double c_k1;
for (int j = n_discrete; j < end; ++j) {
k = j;
c_k1 = distribution_[l].c[k+1];
if (r1 < c_k1) break;
k = j + 1;
c_k = c_k1;
}
double E_l_k = distribution_[l].e_out[k];
double p_l_k = distribution_[l].p[k];
double km_r, km_a;
if (distribution_[l].interpolation == Interpolation::histogram) {
// Histogram interpolation
if (p_l_k > 0.0 && k >= distribution_[l].n_discrete) {
if (p_l_k > 0.0 && k >= n_discrete) {
E_out = E_l_k + (r1 - c_k)/p_l_k;
} else {
E_out = E_l_k;
@ -205,10 +220,12 @@ void KalbachMann::sample(double E_in, double& E_out, double& mu) const
}
// Now interpolate between incident energy bins i and i + 1
if (l == i) {
E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1);
} else {
E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1);
if (k >= n_discrete) {
if (l == i) {
E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1);
} else {
E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1);
}
}
// Sampled correlated angle from Kalbach-Mann parameters