Remove Material type on Fortran side

This commit is contained in:
Paul Romano 2019-01-27 16:10:53 -06:00
parent 63aa3bf106
commit b25feff0e7
9 changed files with 70 additions and 253 deletions

View file

@ -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))

View file

@ -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

View file

@ -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

View file

@ -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 <material>
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()

View file

@ -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();

View file

@ -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

View file

@ -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

View file

@ -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

View file

@ -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