From f407982ff65b7346aa5b2aec58b8c371fefee3d9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 25 Jan 2019 11:17:03 -0600 Subject: [PATCH] Remove Material p0 and element attributes on Fortran side --- include/openmc/material.h | 1 - src/input_xml.F90 | 58 --------------------------------------- src/material.cpp | 17 +++--------- src/material_header.F90 | 29 -------------------- src/physics.cpp | 42 ++++++++++++++-------------- 5 files changed, 25 insertions(+), 122 deletions(-) diff --git a/include/openmc/material.h b/include/openmc/material.h index ee645c287..0443c3480 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -73,7 +73,6 @@ public: double volume_ {-1.0}; //!< Volume in [cm^3] bool fissionable_ {false}; //!< Does this material contain fissionable nuclides bool depletable_ {false}; //!< Is the material depletable? - bool has_isotropic_nuclides_ {false}; std::vector p0_; //!< Indicate which nuclides are to be treated with iso-in-lab scattering // Thermal scattering tables diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 9b8e2f08e..c97642e1f 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -388,22 +388,18 @@ contains subroutine read_materials_xml() bind(C) integer :: i ! loop index for materials integer :: j ! loop index for nuclides - integer :: k ! loop index integer :: n ! number of nuclides integer :: index_nuclide ! index in nuclides integer :: index_element ! index in elements integer :: index_sab ! index in sab_tables logical :: file_exists ! does materials.xml exist? character(20) :: name ! name of nuclide, e.g. U235 - character(3) :: element ! name of element, e.g. Zr character(MAX_WORD_LEN) :: units ! units on density character(MAX_LINE_LEN) :: filename ! absolute path to materials.xml - character(MAX_WORD_LEN), allocatable :: sarray(:) real(8) :: val ! value entered for density real(8) :: temp_dble ! temporary double prec. real logical :: sum_density ! density is sum of nuclide densities type(VectorChar) :: names ! temporary list of nuclide names - type(VectorChar) :: list_iso_lab ! temporary list of isotropic lab scatterers type(VectorReal) :: densities ! temporary list of nuclide densities type(Material), pointer :: mat => null() type(XMLDocument) :: doc @@ -624,19 +620,6 @@ contains end do INDIVIDUAL_NUCLIDES end if - ! ======================================================================= - ! READ AND PARSE element - - if (check_for_node(node_mat, "isotropic")) then - n = node_word_count(node_mat, "isotropic") - allocate(sarray(n)) - call get_node_array(node_mat, "isotropic", sarray) - do j = 1, n - call list_iso_lab % push_back(sarray(j)) - end do - deallocate(sarray) - end if - ! ======================================================================== ! COPY NUCLIDES TO ARRAYS IN MATERIAL @@ -645,7 +628,6 @@ contains mat % n_nuclides = n allocate(mat % names(n)) allocate(mat % nuclide(n)) - allocate(mat % element(n)) allocate(mat % atom_density(n)) ALL_NUCLIDES: do j = 1, mat % n_nuclides @@ -667,51 +649,12 @@ contains mat % nuclide(j) = nuclide_dict % get(name) end if - ! If the corresponding element hasn't been encountered yet and photon - ! transport will be used, we need to add its symbol to the element_dict - if (photon_transport) then - element = name(1:scan(name, '0123456789') - 1) - - ! Make sure photon cross section data is available - if (.not. library_present(LIBRARY_PHOTON, element)) then - call fatal_error("Could not find element " // trim(element) & - // " in cross_sections data file!") - end if - - if (.not. element_dict % has(element)) then - index_element = index_element + 1 - mat % element(j) = index_element - - call element_dict % set(element, index_element) - else - mat % element(j) = element_dict % get(element) - end if - end if - ! Copy name and atom/weight percent mat % names(j) = name mat % atom_density(j) = densities % data(j) end do ALL_NUCLIDES - if (run_CE) then - ! By default, isotropic-in-lab is not used - if (list_iso_lab % size() > 0) then - mat % has_isotropic_nuclides = .true. - allocate(mat % p0(n)) - mat % p0(:) = .false. - - ! Apply isotropic-in-lab treatment to specified nuclides - do j = 1, list_iso_lab % size() - do k = 1, n - if (names % data(k) == list_iso_lab % data(j)) then - mat % p0(k) = .true. - end if - end do - end do - end if - end if - ! Check to make sure either all atom percents or all weight percents are ! given if (.not. (all(mat % atom_density >= ZERO) .or. & @@ -726,7 +669,6 @@ contains ! Clear lists call names % clear() call densities % clear() - call list_iso_lab % clear() ! Add material to dictionary call material_dict % set(mat % id(), i) diff --git a/src/material.cpp b/src/material.cpp index 05e7dc57c..dcd24b806 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -252,7 +252,6 @@ Material::Material(pugi::xml_node node) if (settings::run_CE) { // By default, isotropic-in-lab is not used if (iso_lab.size() > 0) { - has_isotropic_nuclides_ = true; p0_.resize(n); // Apply isotropic-in-lab treatment to specified nuclides @@ -458,14 +457,8 @@ void Material::init_bremsstrahlung() auto n_k = data::ttb_k_grid.size(); auto n_e = data::ttb_e_grid.size(); - // Get pointers to nuclides, elements, densities - int32_t index; - openmc_get_material_index(id_, &index); - int* nuclide_; - double* atom_density_; - int n; - openmc_material_get_densities(index, &nuclide_, &atom_density_, &n); - int* element_ = material_element(index); + // Determine number of elements + int n = element_.size(); for (int particle = 0; particle < 2; ++particle) { // Loop over logic twice, once for electron, once for positron @@ -497,10 +490,8 @@ void Material::init_bremsstrahlung() // fixed in the future. for (int i = 0; i < n; ++i) { // Get pointer to current element - // TODO: off-by-one - const auto& elm = data::elements[element_[i] - 1]; - // TODO: off-by-one - double awr = data::nuclides[nuclide_[i] - 1]->awr_; + const auto& elm = data::elements[element_[i]]; + double awr = data::nuclides[nuclide_[i]]->awr_; // Get atomic density and mass density of nuclide given atom/weight percent double atom_density = (atom_density_[0] > 0.0) ? diff --git a/src/material_header.F90 b/src/material_header.F90 index 84eda36df..8246a082a 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -70,7 +70,6 @@ module material_header character(len=104) :: name = "" ! User-defined name integer :: n_nuclides = 0 ! number of nuclides integer, allocatable :: nuclide(:) ! index in nuclides array - integer, allocatable :: element(:) ! index in elements array real(8) :: density ! total atom density in atom/b-cm real(C_DOUBLE), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm real(8) :: density_gpcc ! total density in g/cm^3 @@ -93,11 +92,6 @@ module material_header ! Does this material contain fissionable nuclides? Is it depletable? logical :: depletable = .false. - - ! enforce isotropic scattering in lab for specific nuclides - logical :: has_isotropic_nuclides = .false. - logical, allocatable :: p0(:) - contains procedure :: id => material_id procedure :: set_id => material_set_id @@ -233,27 +227,4 @@ contains err = 0 end function openmc_extend_materials -!=============================================================================== -! Fortran compatibility -!=============================================================================== - - function material_isotropic(i_material, i_nuc_mat) result(iso) bind(C) - integer(C_INT), value :: i_material - integer(C_INT), value :: i_nuc_mat - logical(C_BOOL) :: iso - - iso = .false. - associate (mat => materials(i_material)) - if (mat % has_isotropic_nuclides) then - iso = mat % p0(i_nuc_mat) - end if - end associate - end function - - function material_element(i_material) result(ptr) bind(C) - integer(C_INT), value :: i_material - type(C_PTR) :: ptr - ptr = C_LOC(materials(i_material) % element(1)) - end function - end module material_header diff --git a/src/physics.cpp b/src/physics.cpp index a00897d4f..770380594 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -483,11 +483,8 @@ int sample_element(Particle* p) double cutoff = prn() * simulation::material_xs.total; // Get pointers to elements, densities - int* nuclide; - double* density; - int n; - openmc_material_get_densities(p->material, &nuclide, &density, &n); - int* element = material_element(p->material); + const auto& mat {model::materials[p->material - 1]}; + int n = mat->nuclide_.size(); int i = 0; double prob = 0.0; @@ -500,9 +497,8 @@ int sample_element(Particle* p) } // Find atom density - // TODO: off-by-one - i_element = element[i] - 1; - double atom_density = density[i]; + i_element = mat->element_[i]; + double atom_density = mat->atom_density_[i]; // Determine microscopic cross section double sigma = atom_density * simulation::micro_photon_xs[i_element].total; @@ -710,21 +706,25 @@ void scatter(Particle* p, int i_nuclide, int i_nuc_mat) p->event = EVENT_SCATTER; // Sample new outgoing angle for isotropic-in-lab scattering - if (material_isotropic(p->material, i_nuc_mat)) { - // Sample isotropic-in-lab outgoing direction - double mu = 2.0*prn() - 1.0; - double phi = 2.0*PI*prn(); - Direction u_new; - u_new.x = mu; - u_new.y = std::sqrt(1.0 - mu*mu)*std::cos(phi); - u_new.z = std::sqrt(1.0 - mu*mu)*std::sin(phi); + // TODO: off-by-one + const auto& mat {model::materials[p->material - 1]}; + if (!mat->p0_.empty()) { + if (mat->p0_[i_nuc_mat]) { + // Sample isotropic-in-lab outgoing direction + double mu = 2.0*prn() - 1.0; + double phi = 2.0*PI*prn(); + Direction u_new; + u_new.x = mu; + u_new.y = std::sqrt(1.0 - mu*mu)*std::cos(phi); + u_new.z = std::sqrt(1.0 - mu*mu)*std::sin(phi); - p->mu = u_old.dot(u_new); + p->mu = u_old.dot(u_new); - // Change direction of particle - p->coord[0].uvw[0] = u_new.x; - p->coord[0].uvw[1] = u_new.y; - p->coord[0].uvw[2] = u_new.z; + // Change direction of particle + p->coord[0].uvw[0] = u_new.x; + p->coord[0].uvw[1] = u_new.y; + p->coord[0].uvw[2] = u_new.z; + } } }