From 2eacd5f72ce8f9cd6ac68d5fb1717dfc7ece134a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 12 Jan 2019 17:04:36 -0600 Subject: [PATCH] Finish TTB implementation on C++ side, still a bug somewhere --- CMakeLists.txt | 1 + include/openmc/bremsstrahlung.h | 4 + include/openmc/material.h | 3 +- include/openmc/physics.h | 2 - src/bremsstrahlung.cpp | 162 ++++++++++++++++++-------------- src/input_xml.F90 | 13 ++- src/material.cpp | 32 ++++--- src/photon.cpp | 6 +- src/physics.cpp | 5 +- 9 files changed, 133 insertions(+), 95 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 7e9e24161..808cc8e2f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -370,6 +370,7 @@ add_library(libopenmc SHARED src/tallies/trigger.F90 src/tallies/trigger_header.F90 src/bank.cpp + src/bremsstrahlung.cpp src/dagmc.cpp src/cell.cpp src/cmfd_solver.cpp diff --git a/include/openmc/bremsstrahlung.h b/include/openmc/bremsstrahlung.h index ac74a2ab4..bd9d012c1 100644 --- a/include/openmc/bremsstrahlung.h +++ b/include/openmc/bremsstrahlung.h @@ -37,6 +37,10 @@ extern xt::xtensor ttb_k_grid; //! reduced energy W/T of emitted phot } // namespace data +//============================================================================== +// Global variables +//============================================================================== + void thick_target_bremsstrahlung(Particle& p, double* E_lost); } // namespace openmc diff --git a/include/openmc/material.h b/include/openmc/material.h index 8df8fa308..d796b08ed 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -47,7 +47,7 @@ public: explicit Material(pugi::xml_node material_node); -private: + //! Initialize bremsstrahlung data void init_bremsstrahlung(); }; @@ -55,6 +55,7 @@ private: // Fortran compatibility //============================================================================== +extern "C" int* material_element(int i_material); extern "C" bool material_isotropic(int i_material, int i_nuc_mat); } // namespace openmc diff --git a/include/openmc/physics.h b/include/openmc/physics.h index 8a3d7c4ef..921494db7 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -82,8 +82,6 @@ void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p); void sample_secondary_photons(Particle* p, int i_nuclide); -extern "C" void thick_target_bremsstrahlung(Particle* p, double* E_lost); - } // namespace openmc #endif // OPENMC_PHYSICS_H diff --git a/src/bremsstrahlung.cpp b/src/bremsstrahlung.cpp index 1ccf1f351..1ba1073a5 100644 --- a/src/bremsstrahlung.cpp +++ b/src/bremsstrahlung.cpp @@ -1,5 +1,13 @@ #include "openmc/bremsstrahlung.h" +#include "openmc/constants.h" +#include "openmc/material.h" +#include "openmc/random_lcg.h" +#include "openmc/search.h" +#include "openmc/settings.h" + +#include "xtensor/xmath.hpp" + namespace openmc { //============================================================================== @@ -18,92 +26,104 @@ std::vector ttb; // Non-member functions //============================================================================== +void thick_target_bremsstrahlung(Particle& p, double* E_lost) +{ + if (p.material == MATERIAL_VOID) return; -// void thick_target_bremsstrahlung(Particle& p, double* E_lost) -// { -// if (p.material == MATERIAL_VOID) return; + // TODO: off-by-one + int photon = static_cast(ParticleType::photon) - 1; + if (p.E < settings::energy_cutoff[photon]) return; -// // TODO: off-by-one -// int photon = static_cast(ParticleType::photon) - 1; -// if (p.E < settings::energy_cutoff[photon]) return; + // Get bremsstrahlung data for this material and particle type + // TODO: off-by-one + BremsstrahlungData* mat; + if (p.type == static_cast(ParticleType::positron)) { + mat = &model::materials[p.material -1]->ttb_->positron; + } else { + mat = &model::materials[p.material -1]->ttb_->electron; + } -// // // Get bremsstrahlung data for this material and particle type -// if (p.type == static_cast(ParticleType::positron)) { -// mat => ttb(p.material) % positron; -// } else { -// mat => ttb(p.material) % electron; -// } + double e = std::log(p.E); + auto n_e = data::ttb_e_grid.size(); -// double e = std::log(p.E); -// auto n_e = data::ttb_e_grid + // Find the lower bounding index of the incident electron energy + int j = lower_bound_index(data::ttb_e_grid.cbegin(), + data::ttb_e_grid.cend(), e); + if (j == n_e - 1) --j; -// // Find the lower bounding index of the incident electron energy -// int j = lower_bound_index(data::ttb_e_grid.cbegin(), data::ttb_e_grid.cend(), e); -// if (j == n_e - 1) --j; + // Get the interpolation bounds + double e_l = data::ttb_e_grid(j); + double e_r = data::ttb_e_grid(j+1); + double y_l = mat->yield(j); + double y_r = mat->yield(j+1); -// // Get the interpolation bounds -// double e_l = data::ttb_e_grid(j); -// double e_r = data::ttb_e_grid(j+1); -// double y_l = mat % yield(j); -// double y_r = mat % yield(j+1); + // Calculate the interpolation weight w_j+1 of the bremsstrahlung energy PDF + // interpolated in log energy, which can be interpreted as the probability + // of index j+1 + double f = (e - e_l)/(e_r - e_l); -// // Calculate the interpolation weight w_j+1 of the bremsstrahlung energy PDF -// // interpolated in log energy, which can be interpreted as the probability -// // of index j+1 -// double f = (e - e_l)/(e_r - e_l); + // Get the photon number yield for the given energy using linear + // interpolation on a log-log scale + double y = std::exp(y_l + (y_r - y_l)*f); -// // Get the photon number yield for the given energy using linear -// // interpolation on a log-log scale -// double y = std::exp(y_l + (y_r - y_l)*f); + // Sample number of secondary bremsstrahlung photons + int n = y + prn(); -// // Sample number of secondary bremsstrahlung photons -// int n = y + prn(); + *E_lost = 0.0; + if (n == 0) return; -// *E_lost = 0.0; -// if (n == 0) return; + // Sample index of the tabulated PDF in the energy grid, j or j+1 + double c_max; + int i_e; + if (prn() <= f || j == 0) { + i_e = j + 1; -// // Sample index of the tabulated PDF in the energy grid, j or j+1 -// double c_max; -// if (prn() <= f || j == 0) { -// int i_e = j + 1; + // Interpolate the maximum value of the CDF at the incoming particle + // energy on a log-log scale + double p_l = mat->pdf(i_e, i_e - 1); + double p_r = mat->pdf(i_e, i_e); + double c_l = mat->cdf(i_e, i_e - 1); + double a = std::log(p_r/p_l)/(e_r - e_l) + 1.0; + c_max = c_l + std::exp(e_l)*p_l/a*(std::exp(a*(e - e_l)) - 1.0); + } else { + i_e = j; -// // Interpolate the maximum value of the CDF at the incoming particle -// // energy on a log-log scale -// double p_l = mat % pdf(i_e-1, i_e); -// double p_r = mat % pdf(i_e, i_e); -// double c_l = mat % cdf(i_e-1, i_e); -// double a = std::log(p_r/p_l)/(e_r - e_l) + 1.0; -// c_max = c_l + std::exp(e_l)*p_l/a*(std::exp(a*(e - e_l)) - 1.0); -// } else { -// int i_e = j; + // Maximum value of the CDF + c_max = mat->cdf(i_e, i_e); + } -// // Maximum value of the CDF -// c_max = mat % cdf(i_e, i_e) -// } + // Sample the energies of the emitted photons + for (int i = 0; i < n; ++i) { + // Generate a random number r and determine the index i for which + // cdf(i) <= r*cdf,max <= cdf(i+1) + double c = prn()*c_max; + int i_w = lower_bound_index(&mat->cdf(i_e, 0), &mat->cdf(i_e, 0) + i_e, c); -// // Sample the energies of the emitted photons -// for (int i = 0; i < n; ++i) { -// // Generate a random number r and determine the index i for which -// // cdf(i) <= r*cdf,max <= cdf(i+1) -// double c = prn()*c_max; -// int i_w = lower_bound_index(mat % cdf(:i_e,i_e), c) + // Sample the photon energy + double w_l = data::ttb_e_grid(i_w); + double w_r = data::ttb_e_grid(i_w + 1); + double p_l = mat->pdf(i_e, i_w); + double p_r = mat->pdf(i_e, i_w + 1); + double c_l = mat->cdf(i_e, i_w); + double a = std::log(p_r/p_l)/(w_r - w_l) + 1.0; + double w = std::exp(w_l)*std::pow(a*(c - c_l)/(std::exp(w_l)*p_l) + 1.0, 1.0/a); -// // Sample the photon energy -// double w_l = data::ttb_e_grid(i_w); -// double w_r = data::ttb_e_grid(i_w+1); -// double p_l = mat % pdf(i_w, i_e); -// double p_r = mat % pdf(i_w+1, i_e); -// double c_l = mat % cdf(i_w, i_e); -// double a = std::log(p_r/p_l)/(w_r - w_l) + 1.0; -// double w = std::exp(w_l)*std::pow(a*(c - c_l)/(std::exp(w_l)*p_l) + 1.0, 1.0/a); + if (w > settings::energy_cutoff[photon]) { + // Create secondary photon + int photon_ = static_cast(ParticleType::photon); + p.create_secondary(p.coord[0].uvw, w, photon_, true); + *E_lost += w; + } + } +} -// if (w > settings::energy_cutoff[photon]) { -// // Create secondary photon -// int photon = static_cast(ParticleType::photon); -// p.create_secondary(p.coord[0].uvw, w, photon, true); -// *E_lost += w; -// } -// } -// } +//============================================================================== +// Fortran compatibility +//============================================================================== + +extern "C" void set_log_ttb_e_grid() +{ + data::ttb_e_grid = xt::log(data::ttb_e_grid); +} } // namespace openmc diff --git a/src/input_xml.F90 b/src/input_xml.F90 index be93a164e..4b6d77494 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2034,6 +2034,15 @@ contains type(SetChar) :: already_read type(SetChar) :: element_already_read + interface + subroutine material_init_bremsstrahlung(ptr) bind(C) + import C_PTR + type(C_PTR), value :: ptr + end subroutine + subroutine set_log_ttb_e_grid() bind(C) + end subroutine + end interface + allocate(nuclides(n_nuclides)) allocate(elements(n_elements)) if (photon_transport .and. electron_treatment == ELECTRON_TTB) then @@ -2127,8 +2136,7 @@ contains ! Generate material bremsstrahlung data for electrons and positrons if (photon_transport .and. electron_treatment == ELECTRON_TTB) then - call bremsstrahlung_init(ttb(i) % electron, i, ELECTRON) - call bremsstrahlung_init(ttb(i) % positron, i, POSITRON) + call material_init_bremsstrahlung(materials(i) % ptr) end if end do @@ -2155,6 +2163,7 @@ contains ! Take logarithm of energies since they are log-log interpolated ttb_e_grid = log(ttb_e_grid) + call set_log_ttb_e_grid() end if ! Set up logarithmic grid for nuclides diff --git a/src/material.cpp b/src/material.cpp index 968094994..acbfc08f5 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -49,8 +49,6 @@ Material::Material(pugi::xml_node node) } } -int* material_element(int index); - void Material::init_bremsstrahlung() { // Create new object @@ -80,12 +78,9 @@ void Material::init_bremsstrahlung() ttb->yield = xt::empty({n_e}); // Allocate temporary arrays - std::array shape {n_e}; xt::xtensor stopping_power_collision({n_e}, 0.0); xt::xtensor stopping_power_radiative({n_e}, 0.0); xt::xtensor dcs({n_e, n_k}, 0.0); - xt::xtensor f({n_e}, 0.0); - xt::xtensor z({n_e}, 0.0); double Z_eq_sq = 0.0; double sum_density = 0.0; @@ -154,31 +149,35 @@ void Material::init_bremsstrahlung() stopping_power_radiative; // Loop over photon energies + xt::xtensor f({n_e}, 0.0); + xt::xtensor z({n_e}, 0.0); for (int i = 0; i < n_e - 1; ++i) { double w = data::ttb_e_grid(i); // Loop over incident particle energies - for (int j = 0; j < n_e; ++j) { + for (int j = i; j < n_e; ++j) { double e = data::ttb_e_grid(j); // Reduced photon energy double k = w / e; // Find the lower bounding index of the reduced photon energy - int i_k = lower_bound_index(data::ttb_k_grid.cbegin(), data::ttb_k_grid.cend(), k); + int i_k = lower_bound_index(data::ttb_k_grid.cbegin(), + data::ttb_k_grid.cend(), k); // Get the interpolation bounds double k_l = data::ttb_k_grid(i_k); double k_r = data::ttb_k_grid(i_k + 1); double x_l = dcs(j, i_k); - double x_r = dcs(j, i_k+1); + double x_r = dcs(j, i_k + 1); // Find the value of the DCS using linear interpolation in reduced // photon energy k - double x = x_l + (k - k_l) * (x_r - x_l) / (k_r - k_l); + double x = x_l + (k - k_l)*(x_r - x_l)/(k_r - k_l); // Ratio of the velocity of the charged particle to the speed of light - double beta = std::sqrt(e*(e + 2.0*MASS_ELECTRON_EV)) / (e + MASS_ELECTRON_EV); + double beta = std::sqrt(e*(e + 2.0*MASS_ELECTRON_EV)) / + (e + MASS_ELECTRON_EV); // Compute the integrand of the PDF f(j) = x / (beta*beta * stopping_power(j) * w); @@ -193,7 +192,7 @@ void Material::init_bremsstrahlung() spline_c(n, &data::ttb_e_grid(i), &f(i), &z(i)); double c = 0.0; - for (int j = 0; j < n_e - 1; ++j) { + for (int j = i; j < n_e - 1; ++j) { c += spline_integrate_c(n, &data::ttb_e_grid(i), &f(i), &z(i), data::ttb_e_grid(j), data::ttb_e_grid(j+1)); @@ -207,7 +206,7 @@ void Material::init_bremsstrahlung() double x_l = std::log(f(i)); double x_r = std::log(f(i+1)); - ttb->pdf(i,i+1) = 0.5*(e_r - e_l)*(std::exp(e_l + x_l) + std::exp(e_r + x_r)); + ttb->pdf(i+1,i) = 0.5*(e_r - e_l)*(std::exp(e_l + x_l) + std::exp(e_r + x_r)); } } @@ -219,7 +218,7 @@ void Material::init_bremsstrahlung() // Loop over photon energies double c = 0.0; - for (int i = 0; i < j - 1; ++i) { + for (int i = 0; i < j; ++i) { // Integrate the CDF from the PDF using the trapezoidal rule in log-log // space double w_l = std::log(data::ttb_e_grid(i)); @@ -236,7 +235,7 @@ void Material::init_bremsstrahlung() } // Use logarithm of number yield since it is log-log interpolated - ttb->yield = xt::where(ttb->yield > 0.0, ttb->yield, -500.0); + ttb->yield = xt::where(ttb->yield > 0.0, xt::log(ttb->yield), -500.0); } } @@ -332,6 +331,11 @@ extern "C" { mat->fissionable = fissionable; } + void material_init_bremsstrahlung(Material* mat) + { + mat->init_bremsstrahlung(); + } + void extend_materials_c(int32_t n) { model::materials.reserve(model::materials.size() + n); diff --git a/src/photon.cpp b/src/photon.cpp index b324f5ff3..3b489e299 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -1,5 +1,6 @@ #include "openmc/photon.h" +#include "openmc/bremsstrahlung.h" #include "openmc/constants.h" #include "openmc/hdf5_interface.h" #include "openmc/particle.h" @@ -24,8 +25,6 @@ namespace openmc { namespace data { xt::xtensor compton_profile_pz; -xt::xtensor ttb_e_grid; -xt::xtensor ttb_k_grid; std::vector elements; @@ -250,7 +249,7 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element) frst, xt::view(s_rad, xt::range(i_grid+1, n_e)))); // Interpolate bremsstrahlung DCS at the cutoff energy and truncate - xt::xtensor dcs = xt::empty({n_e - i_grid, n_k}); + xt::xtensor dcs({n_e - i_grid, n_k}); for (int i = 0; i < n_k; ++i) { y = std::exp(std::log(dcs_(i_grid,i)) + f*(std::log(dcs_(i_grid+1,i)) - std::log(dcs_(i_grid,i)))); @@ -268,6 +267,7 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element) } // Set incident particle energy grid + // TODO: Change to zero when xtensor is updated if (data::ttb_e_grid.size() == 1) { data::ttb_e_grid = electron_energy; } diff --git a/src/physics.cpp b/src/physics.cpp index dbc6252eb..564bc0f5c 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -1,6 +1,7 @@ #include "openmc/physics.h" #include "openmc/bank.h" +#include "openmc/bremsstrahlung.h" #include "openmc/constants.h" #include "openmc/eigenvalue.h" #include "openmc/error.h" @@ -382,7 +383,7 @@ void sample_electron_reaction(Particle* p) if (settings::electron_treatment == ELECTRON_TTB) { double E_lost; - thick_target_bremsstrahlung(p, &E_lost); + thick_target_bremsstrahlung(*p, &E_lost); } p->E = 0.0; @@ -395,7 +396,7 @@ void sample_positron_reaction(Particle* p) if (settings::electron_treatment == ELECTRON_TTB) { double E_lost; - thick_target_bremsstrahlung(p, &E_lost); + thick_target_bremsstrahlung(*p, &E_lost); } // Sample angle isotropically