From b25feff0e71f5da7bc24ab1963b195ced05c9649 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 27 Jan 2019 16:10:53 -0600 Subject: [PATCH] Remove Material type on Fortran side --- openmc/capi/material.py | 6 +- src/api.F90 | 4 +- src/geometry_header.F90 | 2 - src/input_xml.F90 | 15 ----- src/material.cpp | 31 +++++---- src/material_header.F90 | 144 +++------------------------------------- src/tallies/tally.F90 | 115 +++++++++----------------------- src/tracking.F90 | 4 +- src/volume_calc.F90 | 2 +- 9 files changed, 70 insertions(+), 253 deletions(-) diff --git a/openmc/capi/material.py b/openmc/capi/material.py index 4d40be04b..4bf19a1e0 100644 --- a/openmc/capi/material.py +++ b/openmc/capi/material.py @@ -1,5 +1,5 @@ from collections.abc import Mapping -from ctypes import c_int, c_int32, c_double, c_char_p, POINTER +from ctypes import c_int, c_int32, c_double, c_char_p, POINTER, c_size_t from weakref import WeakValueDictionary import numpy as np @@ -48,6 +48,8 @@ _dll.openmc_material_set_id.errcheck = _error_handler _dll.openmc_material_set_volume.argtypes = [c_int32, c_double] _dll.openmc_material_set_volume.restype = c_int _dll.openmc_material_set_volume.errcheck = _error_handler +_dll.n_materials.argtypes = [] +_dll.n_materials.restype = c_size_t class Material(_FortranObjectWithID): @@ -228,7 +230,7 @@ class _MaterialMapping(Mapping): yield Material(index=i + 1).id def __len__(self): - return c_int32.in_dll(_dll, 'n_materials').value + return _dll.n_materials() def __repr__(self): return repr(dict(self)) diff --git a/src/api.F90 b/src/api.F90 index 6cb1d905a..08112d752 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -62,7 +62,6 @@ contains use geometry_header use material_header use photon_header - use sab_header use settings use simulation_header use surface_header @@ -76,6 +75,9 @@ contains subroutine free_memory_source() bind(C) end subroutine + subroutine free_memory_material() bind(C) + end subroutine + subroutine free_memory_mesh() bind(C) end subroutine free_memory_mesh diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 132ec93ab..1ac6ba423 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -5,9 +5,7 @@ module geometry_header use algorithm, only: find use constants, only: K_BOLTZMANN, MATERIAL_VOID use dict_header, only: DictIntInt - use material_header, only: Material, materials use nuclide_header - use sab_header use stl_vector, only: VectorReal use string, only: to_lower diff --git a/src/input_xml.F90 b/src/input_xml.F90 index afbcc72fe..f2a304064 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -17,7 +17,6 @@ module input_xml use message_passing use mgxs_interface use nuclide_header - use multipole_header use output, only: title, header use photon_header use random_lcg, only: prn @@ -386,13 +385,10 @@ contains end subroutine allocate_cells subroutine read_materials_xml() bind(C) - integer :: i ! loop index for materials logical :: file_exists ! does materials.xml exist? character(MAX_LINE_LEN) :: filename ! absolute path to materials.xml - type(Material), pointer :: mat => null() type(XMLDocument) :: doc type(XMLNode) :: root - type(XMLNode), allocatable :: node_mat_list(:) interface function nuclides_size() bind(C) result(n) @@ -423,17 +419,6 @@ contains call read_materials(root % ptr) - ! Get pointer to list of XML - call get_node_list(root, "material", node_mat_list) - - ! Allocate materials array - n_materials = size(node_mat_list) - allocate(materials(n_materials)) - - do i = 1, n_materials - materials(i) % ptr = material_pointer(i - 1) - end do - ! Set total number of nuclides and elements n_nuclides = nuclides_size() n_elements = elements_size() diff --git a/src/material.cpp b/src/material.cpp index ade0f0dc0..8c9ac3f3a 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -1056,6 +1056,7 @@ openmc_material_set_densities(int32_t index, int n, const char** name, const dou // Assign S(a,b) tables mat->init_thermal(); + return 0; } else { set_errmsg("Index in materials array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; @@ -1068,6 +1069,7 @@ openmc_material_set_id(int32_t index, int32_t id) if (index >= 1 && index <= model::materials.size()) { model::materials[index - 1]->id_ = id; model::material_map[id] = index - 1; + return 0; } else { set_errmsg("Index in materials array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; @@ -1092,14 +1094,25 @@ openmc_material_set_volume(int32_t index, double volume) } } +extern "C" int +openmc_extend_materials(int32_t n, int32_t* index_start, int32_t* index_end) +{ + if (index_start) *index_start = model::materials.size(); + if (index_end) *index_end = model::materials.size() + n - 1; + for (int32_t i = 0; i < n; i++) { + model::materials.push_back(new Material()); + } + return 0; +} + //============================================================================== // Fortran compatibility functions //============================================================================== extern "C" { - Material* material_pointer(int32_t indx) {return model::materials[indx];} + size_t n_materials() { return model::materials.size(); } - int32_t material_id(Material* mat) {return mat->id_;} + int32_t material_id(int32_t i_mat) {return model::materials[i_mat - 1]->id_;} int material_nuclide(int32_t i_mat, int idx) { @@ -1126,20 +1139,12 @@ extern "C" { return model::materials[i_mat - 1]->mat_nuclide_index_[i_nuc - 1] + 1; } - void material_calculate_xs_c(Material* mat, const Particle* p) + void material_calculate_xs(const Particle* p) { - mat->calculate_xs(*p); + model::materials[p->material - 1]->calculate_xs(*p); } - void extend_materials_c(int32_t n) - { - model::materials.reserve(model::materials.size() + n); - for (int32_t i = 0; i < n; i++) { - model::materials.push_back(new Material()); - } - } - - void free_memory_material_c() + void free_memory_material() { for (Material *mat : model::materials) {delete mat;} model::materials.clear(); diff --git a/src/material_header.F90 b/src/material_header.F90 index 844e0d731..81baeb861 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -2,41 +2,30 @@ module material_header use, intrinsic :: ISO_C_BINDING - use constants - use dict_header, only: DictIntInt - use error - use nuclide_header use particle_header, only: Particle - use photon_header - use sab_header - use stl_vector, only: VectorReal, VectorInt - use string, only: to_str, to_f_string, to_c_string implicit none private - public :: free_memory_material + public :: material_calculate_xs + public :: material_id public :: material_nuclide public :: material_nuclide_size public :: material_nuclide_index public :: material_atom_density public :: material_density_gpcc - public :: openmc_extend_materials - public :: openmc_material_get_volume - public :: material_pointer interface - function material_pointer(mat_ind) bind(C) result(ptr) - import C_PTR, C_INT32_T - integer(C_INT32_T), intent(in), value :: mat_ind - type(C_PTR) :: ptr - end function material_pointer + function material_id(i_mat) bind(C) result(id) + import C_INT32_T + integer(C_INT32_T), value :: i_mat + integer(C_INT32_T) :: id + end function - function material_id_c(mat_ptr) bind(C, name='material_id') result(id) - import C_PTR, C_INT32_T - type(C_PTR), intent(in), value :: mat_ptr - integer(C_INT32_T) :: id - end function material_id_c + subroutine material_calculate_xs(p) bind(C) + import Particle + type(Particle), intent(in) :: p + end subroutine function material_nuclide(i_mat, idx) bind(C) result(nuc) import C_INT32_T, C_INT @@ -70,117 +59,6 @@ module material_header integer(C_INT32_T), value :: i_mat real(C_DOUBLE) :: density end function - - subroutine extend_materials_c(n) bind(C) - import C_INT32_T - integer(C_INT32_T), intent(in), value :: n - end subroutine extend_materials_c - - function openmc_material_get_volume(index, volume) result(err) bind(C) - import C_INT32_T, C_DOUBLE, C_INT - integer(C_INT32_T), value :: index - real(C_DOUBLE), intent(out) :: volume - integer(C_INT) :: err - end function openmc_material_get_volume end interface -!=============================================================================== -! MATERIAL describes a material by its constituent nuclides -!=============================================================================== - - type, public :: Material - type(C_PTR) :: ptr - contains - procedure :: id => material_id - procedure :: calculate_xs => material_calculate_xs - end type Material - - integer(C_INT32_T), public, bind(C) :: n_materials ! # of materials - - type(Material), public, allocatable, target :: materials(:) - -contains - - function material_id(this) result(id) - class(Material), intent(in) :: this - integer(C_INT32_T) :: id - id = material_id_c(this % ptr) - end function material_id - -!=============================================================================== -! MATERIAL_CALCULATE_XS determines the macroscopic cross sections for the -! material the particle is currently traveling through. -!=============================================================================== - - subroutine material_calculate_xs(this, p) - class(Material), intent(in) :: this - type(Particle), intent(in) :: p - - interface - subroutine material_calculate_xs_c(ptr, p) bind(C) - import C_PTR, Particle - type(C_PTR), value :: ptr - type(Particle), intent(in) :: p - end subroutine - end interface - - call material_calculate_xs_c(this % ptr, p) - end subroutine material_calculate_xs - -!=============================================================================== -! FREE_MEMORY_MATERIAL deallocates global arrays defined in this module -!=============================================================================== - - subroutine free_memory_material() - interface - subroutine free_memory_material_c() bind(C) - end subroutine free_memory_material_c - end interface - call free_memory_material_c() - n_materials = 0 - if (allocated(materials)) deallocate(materials) - end subroutine free_memory_material - -!=============================================================================== -! C API FUNCTIONS -!=============================================================================== - - function openmc_extend_materials(n, index_start, index_end) result(err) bind(C) - ! Extend the materials array by n elements - integer(C_INT32_T), value, intent(in) :: n - integer(C_INT32_T), optional, intent(out) :: index_start - integer(C_INT32_T), optional, intent(out) :: index_end - integer(C_INT) :: err - - integer :: i - type(Material), allocatable :: temp(:) ! temporary materials array - - if (n_materials == 0) then - ! Allocate materials array - allocate(materials(n)) - else - ! Allocate materials array with increased size - allocate(temp(n_materials + n)) - - ! Move original materials to temporary array - temp(1:n_materials) = materials(:) - - ! Move allocation from temporary array - call move_alloc(FROM=temp, TO=materials) - end if - - ! Return indices in materials array - if (present(index_start)) index_start = n_materials + 1 - if (present(index_end)) index_end = n_materials + n - n_materials = n_materials + n - - ! Extend the C++ materials array and get pointers to the C++ objects - call extend_materials_c(n) - do i = n_materials - n, n_materials - materials(i) % ptr = material_pointer(i - 1) - end do - - err = 0 - end function openmc_extend_materials - end module material_header diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 066d715e3..20087eab7 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -1086,13 +1086,11 @@ contains else score = ZERO if (p % material /= MATERIAL_VOID) then - associate (mat => materials(p % material)) do l = 1, material_nuclide_size(p % material) i_nuc = material_nuclide(p % material, l) atom_density_ = material_atom_density(p % material, l) score = score + micro_xs(i_nuc) % reaction(m) * atom_density_ * flux end do - end associate end if end if end if @@ -2005,11 +2003,6 @@ contains integer :: i ! loop index for nuclides in material integer :: i_nuclide ! index in nuclides array real(8) :: atom_density ! atom density of single nuclide in atom/b-cm - type(Material), pointer :: mat - - ! Get pointer to current material. We need this in order to determine what - ! nuclides are in the material - mat => materials(p % material) ! ========================================================================== ! SCORE ALL INDIVIDUAL NUCLIDE REACTION RATES @@ -2203,7 +2196,6 @@ contains real(8) :: filter_weight ! combined weight of all filters real(8) :: atom_density logical :: finished ! found all valid bin combinations - type(Material), pointer :: mat ! A loop over all tallies is necessary because we need to simultaneously ! determine different filter bins for the same tally in order to score to it @@ -2261,9 +2253,6 @@ contains if (i_nuclide > 0) then if (p % material /= MATERIAL_VOID) then - ! Get pointer to current material - mat => materials(p % material) - ! Determine index of nuclide in Material % atom_density array j = material_nuclide_index(p % material, i_nuclide) if (j == 0) cycle NUCLIDE_LOOP @@ -2565,7 +2554,6 @@ contains real(8) :: atom_density ! atom density of single nuclide in atom/b-cm real(8) :: filter_weight ! combined weight of all filters logical :: finished ! found all valid bin combinations - type(Material), pointer :: mat ! Determine track-length estimate of flux flux = p % wgt * distance @@ -2631,9 +2619,6 @@ contains if (i_nuclide > 0) then if (p % material /= MATERIAL_VOID) then - ! Get pointer to current material - mat => materials(p % material) - ! Determine index of nuclide in Material % atom_density array j = material_nuclide_index(p % material, i_nuclide) if (j == 0) cycle NUCLIDE_BIN_LOOP @@ -2717,7 +2702,6 @@ contains ! in atom/b-cm real(8) :: filter_weight ! combined weight of all filters logical :: finished ! found all valid bin combinations - type(Material), pointer :: mat ! Determine collision estimate of flux if (survival_biasing) then @@ -2788,9 +2772,6 @@ contains if (i_nuclide > 0) then if (p % material /= MATERIAL_VOID) then - ! Get pointer to current material - mat => materials(p % material) - ! Determine index of nuclide in Material % atom_density array j = material_nuclide_index(p % material, i_nuclide) if (j == 0) cycle NUCLIDE_BIN_LOOP @@ -3031,7 +3012,7 @@ contains case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, & SCORE_NU_FISSION) - if (materials(p % material) % id() == deriv % diff_material) then + if (material_id(p % material) == deriv % diff_material) then score = score * (flux_deriv + ONE & / material_density_gpcc(p % material)) else @@ -3052,7 +3033,7 @@ contains case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, & SCORE_NU_FISSION) - if (materials(p % material) % id() == deriv % diff_material) then + if (material_id(p % material) == deriv % diff_material) then score = score * (flux_deriv + ONE & / material_density_gpcc(p % material)) else @@ -3094,9 +3075,8 @@ contains case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, & SCORE_NU_FISSION) - if (materials(p % material) % id() == deriv % diff_material & + if (material_id(p % material) == deriv % diff_material & .and. p % event_nuclide == deriv % diff_nuclide) then - associate(mat => materials(p % material)) ! Search for the index of the perturbed nuclide. do l = 1, material_nuclide_size(p % material) if (material_nuclide(p % material, l) == deriv % diff_nuclide) exit @@ -3104,7 +3084,6 @@ contains score = score * (flux_deriv & + ONE / material_atom_density(p % material, l)) - end associate else score = score * flux_deriv end if @@ -3116,7 +3095,7 @@ contains case (ESTIMATOR_COLLISION) scoring_diff_nuclide = & - (materials(p % material) % id() == deriv % diff_material) & + (material_id(p % material) == deriv % diff_material) & .and. (i_nuclide == deriv % diff_nuclide) select case (score_bin) @@ -3126,7 +3105,7 @@ contains case (SCORE_TOTAL) if (i_nuclide == -1 .and. & - materials(p % material) % id() == deriv % diff_material .and. & + material_id(p % material) == deriv % diff_material .and. & material_xs % total /= ZERO) then score = score * (flux_deriv & + micro_xs(deriv % diff_nuclide) % total & @@ -3140,7 +3119,7 @@ contains case (SCORE_SCATTER) if (i_nuclide == -1 .and. & - materials(p % material) % id() == deriv % diff_material .and. & + material_id(p % material) == deriv % diff_material .and. & material_xs % total - material_xs % absorption /= ZERO) then score = score * (flux_deriv & + (micro_xs(deriv % diff_nuclide) % total & @@ -3156,7 +3135,7 @@ contains case (SCORE_ABSORPTION) if (i_nuclide == -1 .and. & - materials(p % material) % id() == deriv % diff_material .and. & + material_id(p % material) == deriv % diff_material .and. & material_xs % absorption /= ZERO) then score = score * (flux_deriv & + micro_xs(deriv % diff_nuclide) % absorption & @@ -3170,7 +3149,7 @@ contains case (SCORE_FISSION) if (i_nuclide == -1 .and. & - materials(p % material) % id() == deriv % diff_material .and. & + material_id(p % material) == deriv % diff_material .and. & material_xs % fission /= ZERO) then score = score * (flux_deriv & + micro_xs(deriv % diff_nuclide) % fission & @@ -3184,7 +3163,7 @@ contains case (SCORE_NU_FISSION) if (i_nuclide == -1 .and. & - materials(p % material) % id() == deriv % diff_material .and. & + material_id(p % material) == deriv % diff_material .and. & material_xs % nu_fission /= ZERO) then score = score * (flux_deriv & + micro_xs(deriv % diff_nuclide) % nu_fission & @@ -3230,9 +3209,8 @@ contains score = score * flux_deriv case (SCORE_TOTAL) - if (materials(p % material) % id() == deriv % diff_material .and. & + if (material_id(p % material) == deriv % diff_material .and. & micro_xs(p % event_nuclide) % total > ZERO) then - associate(mat => materials(p % material)) ! Search for the index of the perturbed nuclide. do l = 1, material_nuclide_size(p % material) if (material_nuclide(p % material, l) == p % event_nuclide) exit @@ -3249,16 +3227,14 @@ contains score = score * (flux_deriv & + (dsig_s + dsig_a) * material_atom_density(p % material, l) & / material_xs % total) - end associate else score = score * flux_deriv end if case (SCORE_SCATTER) - if (materials(p % material) % id() == deriv % diff_material .and. & + if (material_id(p % material) == deriv % diff_material .and. & (micro_xs(p % event_nuclide) % total & - micro_xs(p % event_nuclide) % absorption) > ZERO) then - associate(mat => materials(p % material)) ! Search for the index of the perturbed nuclide. do l = 1, material_nuclide_size(p % material) if (material_nuclide(p % material, l) == p % event_nuclide) exit @@ -3273,15 +3249,13 @@ contains end associate score = score * (flux_deriv + dsig_s * material_atom_density(p % material, l) / & (material_xs % total - material_xs % absorption)) - end associate else score = score * flux_deriv end if case (SCORE_ABSORPTION) - if (materials(p % material) % id() == deriv % diff_material .and. & + if (material_id(p % material) == deriv % diff_material .and. & micro_xs(p % event_nuclide) % absorption > ZERO) then - associate(mat => materials(p % material)) ! Search for the index of the perturbed nuclide. do l = 1, material_nuclide_size(p % material) if (material_nuclide(p % material, l) == p % event_nuclide) exit @@ -3296,15 +3270,13 @@ contains end associate score = score * (flux_deriv + dsig_a * material_atom_density(p % material, l) & / material_xs % absorption) - end associate else score = score * flux_deriv end if case (SCORE_FISSION) - if (materials(p % material) % id() == deriv % diff_material .and. & + if (material_id(p % material) == deriv % diff_material .and. & micro_xs(p % event_nuclide) % fission > ZERO) then - associate(mat => materials(p % material)) ! Search for the index of the perturbed nuclide. do l = 1, material_nuclide_size(p % material) if (material_nuclide(p % material, l) == p % event_nuclide) exit @@ -3319,15 +3291,13 @@ contains end associate score = score * (flux_deriv & + dsig_f * material_atom_density(p % material, l) / material_xs % fission) - end associate else score = score * flux_deriv end if case (SCORE_NU_FISSION) - if (materials(p % material) % id() == deriv % diff_material .and. & + if (material_id(p % material) == deriv % diff_material .and. & micro_xs(p % event_nuclide) % nu_fission > ZERO) then - associate(mat => materials(p % material)) ! Search for the index of the perturbed nuclide. do l = 1, material_nuclide_size(p % material) if (material_nuclide(p % material, l) == p % event_nuclide) exit @@ -3344,7 +3314,6 @@ contains + dsig_f * material_atom_density(p % material, l) / material_xs % nu_fission& * micro_xs(p % event_nuclide) % nu_fission & / micro_xs(p % event_nuclide) % fission) - end associate else score = score * flux_deriv end if @@ -3363,10 +3332,9 @@ contains case (SCORE_TOTAL) if (i_nuclide == -1 .and. & - materials(p % material) % id() == deriv % diff_material .and. & + material_id(p % material) == deriv % diff_material .and. & material_xs % total > ZERO) then cum_dsig = ZERO - associate(mat => materials(p % material)) do l = 1, material_nuclide_size(p % material) i_nuc = material_nuclide(p % material, l) associate (nuc => nuclides(i_nuc)) @@ -3379,10 +3347,9 @@ contains end if end associate end do - end associate score = score * (flux_deriv & + cum_dsig / material_xs % total) - else if (materials(p % material) % id() == deriv % diff_material & + else if (material_id(p % material) == deriv % diff_material & .and. material_xs % total > ZERO) then dsig_s = ZERO dsig_a = ZERO @@ -3400,10 +3367,9 @@ contains case (SCORE_SCATTER) if (i_nuclide == -1 .and. & - materials(p % material) % id() == deriv % diff_material .and. & + material_id(p % material) == deriv % diff_material .and. & (material_xs % total - material_xs % absorption) > ZERO) then cum_dsig = ZERO - associate(mat => materials(p % material)) do l = 1, material_nuclide_size(p % material) i_nuc = material_nuclide(p % material, l) associate (nuc => nuclides(i_nuc)) @@ -3416,10 +3382,9 @@ contains end if end associate end do - end associate score = score * (flux_deriv + cum_dsig & / (material_xs % total - material_xs % absorption)) - else if ( materials(p % material) % id() == deriv % diff_material & + else if ( material_id(p % material) == deriv % diff_material & .and. (material_xs % total - material_xs % absorption) > ZERO)& then dsig_s = ZERO @@ -3438,10 +3403,9 @@ contains case (SCORE_ABSORPTION) if (i_nuclide == -1 .and. & - materials(p % material) % id() == deriv % diff_material .and. & + material_id(p % material) == deriv % diff_material .and. & material_xs % absorption > ZERO) then cum_dsig = ZERO - associate(mat => materials(p % material)) do l = 1, material_nuclide_size(p % material) i_nuc = material_nuclide(p % material, l) associate (nuc => nuclides(i_nuc)) @@ -3453,10 +3417,9 @@ contains end if end associate end do - end associate score = score * (flux_deriv & + cum_dsig / material_xs % absorption) - else if (materials(p % material) % id() == deriv % diff_material & + else if (material_id(p % material) == deriv % diff_material & .and. material_xs % absorption > ZERO) then dsig_a = ZERO associate (nuc => nuclides(i_nuclide)) @@ -3473,10 +3436,9 @@ contains case (SCORE_FISSION) if (i_nuclide == -1 .and. & - materials(p % material) % id() == deriv % diff_material .and. & + material_id(p % material) == deriv % diff_material .and. & material_xs % fission > ZERO) then cum_dsig = ZERO - associate(mat => materials(p % material)) do l = 1, material_nuclide_size(p % material) i_nuc = material_nuclide(p % material, l) associate (nuc => nuclides(i_nuc)) @@ -3488,10 +3450,9 @@ contains end if end associate end do - end associate score = score * (flux_deriv & + cum_dsig / material_xs % fission) - else if (materials(p % material) % id() == deriv % diff_material & + else if (material_id(p % material) == deriv % diff_material & .and. material_xs % fission > ZERO) then dsig_f = ZERO associate (nuc => nuclides(i_nuclide)) @@ -3508,10 +3469,9 @@ contains case (SCORE_NU_FISSION) if (i_nuclide == -1 .and. & - materials(p % material) % id() == deriv % diff_material .and. & + material_id(p % material) == deriv % diff_material .and. & material_xs % nu_fission > ZERO) then cum_dsig = ZERO - associate(mat => materials(p % material)) do l = 1, material_nuclide_size(p % material) i_nuc = material_nuclide(p % material, l) associate (nuc => nuclides(i_nuc)) @@ -3525,10 +3485,9 @@ contains end if end associate end do - end associate score = score * (flux_deriv & + cum_dsig / material_xs % nu_fission) - else if (materials(p % material) % id() == deriv % diff_material & + else if (material_id(p % material) == deriv % diff_material & .and. material_xs % nu_fission > ZERO) then dsig_f = ZERO associate (nuc => nuclides(i_nuclide)) @@ -3576,30 +3535,25 @@ contains select case (deriv % variable) case (DIFF_DENSITY) - associate (mat => materials(p % material)) - if (mat % id() == deriv % diff_material) then + if (material_id(p % material) == deriv % diff_material) then ! phi is proportional to e^(-Sigma_tot * dist) ! (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist ! (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist deriv % flux_deriv = deriv % flux_deriv & - - distance * material_xs % total / material_density_gpcc(p % material) + - distance * material_xs % total / material_density_gpcc(p % material) end if - end associate case (DIFF_NUCLIDE_DENSITY) - associate (mat => materials(p % material)) - if (mat % id() == deriv % diff_material) then + if (material_id(p % material) == deriv % diff_material) then ! phi is proportional to e^(-Sigma_tot * dist) ! (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist ! (1 / phi) * (d_phi / d_N) = - sigma_tot * dist deriv % flux_deriv = deriv % flux_deriv & - distance * micro_xs(deriv % diff_nuclide) % total end if - end associate case (DIFF_TEMPERATURE) - associate (mat => materials(p % material)) - if (mat % id() == deriv % diff_material) then + if (material_id(p % material) == deriv % diff_material) then do l=1, material_nuclide_size(p % material) associate (nuc => nuclides(material_nuclide(p % material, l))) if (multipole_in_range(nuc % ptr, p % E)) then @@ -3614,7 +3568,6 @@ contains end associate end do end if - end associate end select end associate end do @@ -3650,19 +3603,16 @@ contains select case (deriv % variable) case (DIFF_DENSITY) - associate (mat => materials(p % material)) - if (mat % id() == deriv % diff_material) then + if (material_id(p % material) == deriv % diff_material) then ! phi is proportional to Sigma_s ! (1 / phi) * (d_phi / d_rho) = (d_Sigma_s / d_rho) / Sigma_s ! (1 / phi) * (d_phi / d_rho) = 1 / rho deriv % flux_deriv = deriv % flux_deriv & + ONE / material_density_gpcc(p % material) end if - end associate case (DIFF_NUCLIDE_DENSITY) - associate (mat => materials(p % material)) - if (mat % id() == deriv % diff_material & + if (material_id(p % material) == deriv % diff_material & .and. p % event_nuclide == deriv % diff_nuclide) then ! Find the index in this material for the diff_nuclide. do j = 1, material_nuclide_size(p % material) @@ -3679,11 +3629,9 @@ contains deriv % flux_deriv = deriv % flux_deriv & + ONE / material_atom_density(p % material, j) end if - end associate case (DIFF_TEMPERATURE) - associate (mat => materials(p % material)) - if (mat % id() == deriv % diff_material) then + if (material_id(p % material) == deriv % diff_material) then do l=1, material_nuclide_size(p % material) i_nuc = material_nuclide(p % material, l) associate (nuc => nuclides(i_nuc)) @@ -3708,7 +3656,6 @@ contains end associate end do end if - end associate end select end associate end do diff --git a/src/tracking.F90 b/src/tracking.F90 index bf8332852..43df27715 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -11,7 +11,7 @@ module tracking use geometry, only: next_cell #endif - use material_header, only: materials, Material + use material_header, only: material_calculate_xs use message_passing use mgxs_interface use nuclide_header @@ -130,7 +130,7 @@ contains ! If the material is the same as the last material and the ! temperature hasn't changed, we don't need to lookup cross ! sections again. - call materials(p % material) % calculate_xs(p) + call material_calculate_xs(p) end if else ! Get the MG data diff --git a/src/volume_calc.F90 b/src/volume_calc.F90 index ff079c658..55c872cb3 100644 --- a/src/volume_calc.F90 +++ b/src/volume_calc.F90 @@ -215,7 +215,7 @@ contains i_material = p % material if (i_material /= MATERIAL_VOID) then do i_domain = 1, size(this % domain_id) - if (materials(i_material) % id() == this % domain_id(i_domain)) then + if (material_id(i_material) == this % domain_id(i_domain)) then call check_hit(i_domain, i_material, indices, hits, n_mat) end if end do