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

Almost complete mesh conversion

Browse source
This commit is contained in:
Paul Romano 2018-08-30 14:22:55 -05:00
parent 21a79eefd9
commit 04e48224e6
23 changed files with 380 additions and 446 deletions
Download patch file Download diff file Expand all files Collapse all files

1
CMakeLists.txt
View file

@ -385,6 +385,7 @@ add_library(libopenmc SHARED
src/distribution_energy.cpp
src/distribution_multi.cpp
src/distribution_spatial.cpp
src/eigenvalue.cpp
src/endf.cpp
src/initialize.cpp
src/finalize.cpp

4
include/openmc/eigenvalue.h
View file

@ -2,6 +2,9 @@
#define OPENMC_EIGENVALUE_H
#include <cstdint> // for int64_t
#include <vector>
#include "xtensor/xtensor.hpp"
#include "openmc/particle.h"
@ -15,6 +18,7 @@ extern std::vector<double> entropy; //!< Shannon entropy at each generation
extern xt::xtensor<double, 1> source_frac; //!< Source fraction for UFS
extern "C" int64_t n_bank;
#pragma omp threadprivate(n_bank)
//==============================================================================
// Non-member functions

8
include/openmc/hdf5_interface.h
View file

@ -317,6 +317,14 @@ write_attribute(hid_t obj_id, const char* name, const std::array<T, N>& buffer)
write_attr(obj_id, 1, dims, name, H5TypeMap<T>::type_id, buffer.data());
}
template<typename T> inline void
write_attribute(hid_t obj_id, const char* name, const std::vector<T>& buffer)
{
hsize_t dims[] {buffer.size()};
write_attr(obj_id, 1, dims, name, H5TypeMap<T>::type_id, buffer.data());
}
//==============================================================================
// Templates/overloads for write_dataset
//==============================================================================

2
include/openmc/particle.h
View file

@ -161,7 +161,7 @@ extern "C" {
{mark_as_lost(message.str());}
//! create a particle restart HDF5 file
void write_restart();
void write_restart() const;
};

7
src/cmfd_data.F90
View file

@ -80,7 +80,7 @@ contains
integer :: i_mesh ! flattend index for mesh
logical :: energy_filters! energy filters present
real(8) :: flux ! temp variable for flux
type(RegularMesh), pointer :: m ! pointer for mesh object
type(RegularMesh) :: m ! pointer for mesh object
! Extract spatial and energy indices from object
nx = cmfd % indices(1)
@ -99,7 +99,7 @@ contains
select type(filt => filters(i_filter_mesh) % obj)
type is (MeshFilter)
m => meshes(filt % mesh)
m = meshes(filt % mesh)
end select
! Set mesh widths
@ -354,9 +354,6 @@ contains
! Normalize openmc source distribution
cmfd % openmc_src = cmfd % openmc_src/sum(cmfd % openmc_src)*cmfd%norm
! Nullify all pointers
if (associated(m)) nullify(m)
end subroutine compute_xs
!===============================================================================

3
src/cmfd_header.F90
View file

@ -95,7 +95,8 @@ module cmfd_header
! Main object
type(cmfd_type), public :: cmfd
type(RegularMesh), public, pointer :: cmfd_mesh => null()
integer, public :: index_cmfd_mesh
type(RegularMesh), public :: cmfd_mesh
! Pointers for different tallies
type(TallyContainer), public, pointer :: cmfd_tallies(:) => null()

110
src/cmfd_input.F90
View file

@ -3,7 +3,7 @@ module cmfd_input
use, intrinsic :: ISO_C_BINDING
use cmfd_header
use mesh_header, only: mesh_dict
use mesh_header
use mgxs_interface, only: energy_bins, num_energy_groups
use tally
use tally_header
@ -241,7 +241,7 @@ contains
use constants, only: MAX_LINE_LEN
use error, only: fatal_error, warning
use mesh_header, only: RegularMesh, openmc_extend_meshes
use mesh_header
use string
use tally, only: openmc_tally_allocate
use tally_header, only: openmc_extend_tallies
@ -266,112 +266,22 @@ contains
integer :: iarray3(3) ! temp integer array
real(8) :: rarray3(3) ! temp double array
real(C_DOUBLE), allocatable :: energies(:)
type(RegularMesh), pointer :: m
type(XMLNode) :: node_mesh
err = openmc_extend_meshes(1, i_start)
! Read CMFD mesh
call read_meshes(root % ptr)
! Allocate mesh
cmfd_mesh => meshes(i_start)
m => meshes(i_start)
! Get index of cmfd mesh and set ID
i_start = n_meshes() - 1
err = openmc_mesh_set_id(i_start, i_start)
! Set mesh id
m % id = i_start
! Set mesh type to rectangular
m % type = MESH_REGULAR
! Save reference to CMFD mesh
index_cmfd_mesh = i_start
cmfd_mesh = meshes(i_start)
! Get pointer to mesh XML node
node_mesh = root % child("mesh")
! Determine number of dimensions for mesh
n = node_word_count(node_mesh, "dimension")
if (n /= 2 .and. n /= 3) then
call fatal_error("Mesh must be two or three dimensions.")
end if
m % n_dimension = n
! Allocate attribute arrays
allocate(m % dimension(n))
allocate(m % lower_left(n))
allocate(m % width(n))
allocate(m % upper_right(n))
! Check that dimensions are all greater than zero
call get_node_array(node_mesh, "dimension", iarray3(1:n))
if (any(iarray3(1:n) <= 0)) then
call fatal_error("All entries on the <dimension> element for a tally mesh&
& must be positive.")
end if
! Read dimensions in each direction
m % dimension = iarray3(1:n)
! Read mesh lower-left corner location
if (m % n_dimension /= node_word_count(node_mesh, "lower_left")) then
call fatal_error("Number of entries on <lower_left> must be the same as &
&the number of entries on <dimension>.")
end if
call get_node_array(node_mesh, "lower_left", m % lower_left)
! Make sure both upper-right or width were specified
if (check_for_node(node_mesh, "upper_right") .and. &
check_for_node(node_mesh, "width")) then
call fatal_error("Cannot specify both <upper_right> and <width> on a &
&tally mesh.")
end if
! Make sure either upper-right or width was specified
if (.not.check_for_node(node_mesh, "upper_right") .and. &
.not.check_for_node(node_mesh, "width")) then
call fatal_error("Must specify either <upper_right> and <width> on a &
&tally mesh.")
end if
if (check_for_node(node_mesh, "width")) then
! Check to ensure width has same dimensions
if (node_word_count(node_mesh, "width") /= &
node_word_count(node_mesh, "lower_left")) then
call fatal_error("Number of entries on <width> must be the same as the &
&number of entries on <lower_left>.")
end if
! Check for negative widths
call get_node_array(node_mesh, "width", rarray3(1:n))
if (any(rarray3(1:n) < ZERO)) then
call fatal_error("Cannot have a negative <width> on a tally mesh.")
end if
! Set width and upper right coordinate
m % width = rarray3(1:n)
m % upper_right = m % lower_left + m % dimension * m % width
elseif (check_for_node(node_mesh, "upper_right")) then
! Check to ensure width has same dimensions
if (node_word_count(node_mesh, "upper_right") /= &
node_word_count(node_mesh, "lower_left")) then
call fatal_error("Number of entries on <upper_right> must be the same &
&as the number of entries on <lower_left>.")
end if
! Check that upper-right is above lower-left
call get_node_array(node_mesh, "upper_right", rarray3(1:n))
if (any(rarray3(1:n) < m % lower_left)) then
call fatal_error("The <upper_right> coordinates must be greater than &
&the <lower_left> coordinates on a tally mesh.")
end if
! Set upper right coordinate and width
m % upper_right = rarray3(1:n)
m % width = (m % upper_right - m % lower_left) / real(m % dimension, 8)
end if
! Set volume fraction
m % volume_frac = ONE/real(product(m % dimension),8)
! Add mesh to dictionary
call mesh_dict % set(m % id, i_start)
! Determine number of filters
energy_filters = check_for_node(node_mesh, "energy")
n = merge(5, 3, energy_filters)

1
src/eigenvalue.F90
View file

@ -7,7 +7,6 @@ module eigenvalue
use error, only: fatal_error, warning
use math, only: t_percentile
use mesh, only: count_bank_sites
use mesh_header, only: RegularMesh, meshes
use message_passing
use random_lcg, only: prn, set_particle_seed, advance_prn_seed
use settings

5
src/eigenvalue.cpp
View file

@ -147,5 +147,10 @@ extern "C" double entropy_c(int i)
return entropy.at(i - 1);
}
extern "C" double entropy_clear()
{
entropy.clear();
}
} // namespace openmc

12
src/input_xml.F90
View file

@ -4,7 +4,7 @@ module input_xml
use algorithm, only: find
use cmfd_input, only: configure_cmfd
use cmfd_header, only: cmfd_mesh
use cmfd_header, only: index_cmfd_mesh
use constants
use dict_header, only: DictIntInt, DictCharInt, DictEntryCI
use endf, only: reaction_name
@ -86,11 +86,6 @@ module input_xml
type(C_PTR), value :: node_ptr
end subroutine read_materials
subroutine read_meshes(node_ptr) bind(C)
import C_PTR
type(C_PTR), value :: node_ptr
end subroutine
function find_root_universe() bind(C) result(root)
import C_INT32_T
integer(C_INT32_T) :: root
@ -232,7 +227,6 @@ contains
if (run_mode == MODE_EIGENVALUE) then
! Preallocate space for keff and entropy by generation
call k_generation % reserve(n_max_batches*gen_per_batch)
call entropy % reserve(n_max_batches*gen_per_batch)
end if
! Particle tracks
@ -1222,7 +1216,7 @@ contains
! READ MESH DATA
! Check for user meshes and allocate
call read_meshes()
call read_meshes(root % ptr)
! We only need the mesh info for plotting
if (run_mode == MODE_PLOTTING) then
@ -2396,7 +2390,7 @@ contains
&meshlines on plot " // trim(to_str(pl % id)))
end if
pl % meshlines_mesh => cmfd_mesh
pl % index_meshlines_mesh = index_cmfd_mesh
case ('entropy')

6
src/mesh.cpp
View file

@ -532,10 +532,10 @@ xt::xarray<double> RegularMesh::count_sites(int64_t n, const Bank* bank,
int n_energy, const double* energies, bool* outside)
{
// Determine shape of array for counts
int m = xt::prod(shape_)();
std::size_t m = xt::prod(shape_)();
std::vector<std::size_t> shape;
if (n_energy > 0) {
shape = {m, n_energy};
shape = {m, static_cast<std::size_t>(n_energy)};
} else {
shape = {m};
}
@ -757,6 +757,8 @@ void read_meshes(pugi::xml_node* root)
//==============================================================================
extern "C" {
int n_meshes() { return meshes.size(); }
RegularMesh* mesh_ptr(int i) { return meshes[i-1].get(); }
int32_t mesh_id(RegularMesh* m) { return m->id_; }

42
src/mesh_header.F90
View file

@ -2,35 +2,20 @@ module mesh_header
use, intrinsic :: ISO_C_BINDING
use constants
use dict_header, only: DictIntInt
use error
use hdf5_interface
use string, only: to_str, to_lower
use xml_interface
implicit none
private
public :: openmc_extend_meshes
public :: openmc_get_mesh_index
public :: openmc_mesh_get_id
public :: openmc_mesh_get_dimension
public :: openmc_mesh_get_params
public :: openmc_mesh_set_id
public :: openmc_mesh_set_dimension
public :: openmc_mesh_set_params
!===============================================================================
! STRUCTUREDMESH represents a tessellation of n-dimensional Euclidean space by
! congruent squares or cubes
!===============================================================================
type, public :: RegularMesh
type :: RegularMesh
type(C_PTR) :: ptr
contains
procedure :: id => regular_id
procedure :: volume_frac => regular_volume_frac
procedure :: n_dimension => regular_n_dimension
procedure :: dimension => regular_dimension
procedure :: lower_left => regular_lower_left
procedure :: upper_right => regular_upper_right
procedure :: width => regular_width
@ -41,11 +26,6 @@ module mesh_header
procedure :: get_indices_from_bin => regular_get_indices_from_bin
end type RegularMesh
integer(C_INT32_T), public, bind(C) :: n_meshes = 0 ! # of structured meshes
! Dictionary that maps user IDs to indices in 'meshes'
type(DictIntInt), public :: mesh_dict
interface
function openmc_extend_meshes(n, index_start, index_end) result(err) bind(C)
import C_INT32_T, C_INT
@ -158,21 +138,21 @@ module mesh_header
real(C_DOUBLE) :: w
end function
function mesh_get_bin(ptr, xyz) result(bin) bind(C)
pure function mesh_get_bin(ptr, xyz) result(bin) bind(C)
import C_PTR, C_DOUBLE, C_INT
type(C_PTR), value :: ptr
real(C_DOUBLE), intent(in) :: xyz(*)
integer(C_INT) :: bin
end function
function mesh_get_bin_from_indices(ptr, ijk) result(bin) bind(C)
pure function mesh_get_bin_from_indices(ptr, ijk) result(bin) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: ptr
integer(C_INT), intent(in) :: ijk(*)
integer(C_INT) :: bin
end function
subroutine mesh_get_indices(ptr, xyz, ijk, in_mesh) bind(C)
pure subroutine mesh_get_indices(ptr, xyz, ijk, in_mesh) bind(C)
import C_PTR, C_DOUBLE, C_INT, C_BOOL
type(C_PTR), value :: ptr
real(C_DOUBLE), intent(in) :: xyz(*)
@ -180,7 +160,7 @@ module mesh_header
logical(C_BOOL), intent(out) :: in_mesh
end subroutine
subroutine mesh_get_indices_from_bin(ptr, bin, ijk) bind(C)
pure subroutine mesh_get_indices_from_bin(ptr, bin, ijk) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: ptr
integer(C_INT), value :: bin
@ -192,6 +172,16 @@ module mesh_header
integer(C_INT), value :: i
type(C_PTR) :: ptr
end function
subroutine read_meshes(node_ptr) bind(C)
import C_PTR
type(C_PTR), value :: node_ptr
end subroutine
function n_meshes() result(n) bind(C)
import C_INT
integer(C_INT) :: n
end function
end interface
contains

2
src/particle.cpp
View file

@ -141,7 +141,7 @@ Particle::mark_as_lost(const char* message)
}
void
Particle::write_restart()
Particle::write_restart() const
{
// Dont write another restart file if in particle restart mode
if (settings::run_mode == RUN_MODE_PARTICLE) return;

1
src/physics.F90
View file

@ -6,7 +6,6 @@ module physics
use error, only: fatal_error, warning, write_message
use material_header, only: Material, materials
use math
use mesh_header, only: meshes
use message_passing
use nuclide_header
use particle_header

1
src/physics_mg.F90
View file

@ -7,7 +7,6 @@ module physics_mg
use error, only: fatal_error, warning, write_message
use material_header, only: Material, materials
use math, only: rotate_angle
use mesh_header, only: meshes
use mgxs_interface
use message_passing
use nuclide_header, only: material_xs

6
src/plot.F90
View file

@ -9,7 +9,7 @@ module plot
use hdf5_interface
use output, only: time_stamp
use material_header, only: materials
use mesh_header, only: meshes
use mesh_header, only: meshes, RegularMesh
use particle_header
use plot_header
use progress_header, only: ProgressBar
@ -242,8 +242,8 @@ contains
width = xyz_ur_plot - xyz_ll_plot
m = meshes(pl % index_meshlines_mesh)
call m % get_indices(xyz_ll_plot, ijk_ll(:m % n_dimension), in_mesh)
call m % get_indices(xyz_ur_plot, ijk_ur(:m % n_dimension), in_mesh)
call m % get_indices(xyz_ll_plot, ijk_ll(:m % n_dimension()), in_mesh)
call m % get_indices(xyz_ur_plot, ijk_ur(:m % n_dimension()), in_mesh)
! sweep through all meshbins on this plane and draw borders
do i = ijk_ll(outer), ijk_ur(outer)

1
src/settings.cpp
View file

@ -16,6 +16,7 @@
#include "openmc/error.h"
#include "openmc/file_utils.h"
#include "openmc/mesh.h"
#include "openmc/output.h"
#include "openmc/random_lcg.h"
#include "openmc/source.h"
#include "openmc/string_utils.h"

5
src/simulation.F90
View file

@ -10,8 +10,7 @@ module simulation
use cmfd_execute, only: cmfd_init_batch, cmfd_tally_init, execute_cmfd
use cmfd_header, only: cmfd_on
use constants, only: ZERO
use eigenvalue, only: calculate_average_keff, &
calculate_generation_keff, shannon_entropy, &
use eigenvalue, only: calculate_average_keff, calculate_generation_keff, &
synchronize_bank, keff_generation, k_sum
#ifdef _OPENMP
use eigenvalue, only: join_bank_from_threads
@ -479,7 +478,7 @@ contains
! will potentially populate k_generation and entropy)
current_batch = 0
call k_generation % clear()
call entropy % clear()
call entropy_clear()
need_depletion_rx = .false.
! If this is a restart run, load the state point data and binary source

10
src/simulation_header.F90
View file

@ -64,6 +64,12 @@ module simulation_header
!$omp threadprivate(trace, thread_id, current_work)
interface
subroutine entropy_clear() bind(C)
end subroutine
end interface
contains
!===============================================================================
@ -80,12 +86,12 @@ contains
!===============================================================================
subroutine free_memory_simulation()
if (allocated(work_index)) deallocate(work_index)
call k_generation % clear()
call k_generation % shrink_to_fit()
call entropy % clear()
call entropy % shrink_to_fit()
call entropy_clear()
end subroutine free_memory_simulation
end module simulation_header

2
src/state_point.F90
View file

@ -173,7 +173,7 @@ contains
call write_dataset(file_id, "generations_per_batch", gen_per_batch)
k = k_generation % size()
call write_dataset(file_id, "k_generation", k_generation % data(1:k))
call entropy_to_hdf5()
call entropy_to_hdf5(file_id)
call write_dataset(file_id, "k_col_abs", k_col_abs)
call write_dataset(file_id, "k_col_tra", k_col_tra)
call write_dataset(file_id, "k_abs_tra", k_abs_tra)

277
src/tallies/tally_filter_mesh.F90
View file

@ -92,155 +92,154 @@ contains
real(8) :: distance ! distance traveled in mesh cell
logical :: start_in_mesh ! starting coordinates inside mesh?
logical :: end_in_mesh ! ending coordinates inside mesh?
type(RegularMesh), pointer :: m
type(RegularMesh) :: m
! Get a pointer to the mesh.
m => meshes(this % mesh)
n = m % n_dimension
m = meshes(this % mesh)
if (estimator /= ESTIMATOR_TRACKLENGTH) then
! If this is an analog or collision tally, then there can only be one
! valid mesh bin.
call m % get_bin(p % coord(1) % xyz, bin)
if (bin /= NO_BIN_FOUND) then
if (bin >= 0) then
call match % bins % push_back(bin)
call match % weights % push_back(ONE)
end if
return
end if
! A track can span multiple mesh bins so we need to handle a lot of
! intersection logic for tracklength tallies.
! ! A track can span multiple mesh bins so we need to handle a lot of
! ! intersection logic for tracklength tallies.
! ========================================================================
! Determine if the track intersects the tally mesh.
! ! ========================================================================
! ! Determine if the track intersects the tally mesh.
! Copy the starting and ending coordinates of the particle. Offset these
! just a bit for the purposes of determining if there was an intersection
! in case the mesh surfaces coincide with lattice/geometric surfaces which
! might produce finite-precision errors.
xyz0 = p % last_xyz + TINY_BIT * p % coord(1) % uvw
xyz1 = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw
! ! Copy the starting and ending coordinates of the particle. Offset these
! ! just a bit for the purposes of determining if there was an intersection
! ! in case the mesh surfaces coincide with lattice/geometric surfaces which
! ! might produce finite-precision errors.
! xyz0 = p % last_xyz + TINY_BIT * p % coord(1) % uvw
! xyz1 = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw
! Determine indices for starting and ending location.
call m % get_indices(xyz0, ijk0(:n), start_in_mesh)
call m % get_indices(xyz1, ijk1(:n), end_in_mesh)
! ! Determine indices for starting and ending location.
! call m % get_indices(xyz0, ijk0(:n), start_in_mesh)
! call m % get_indices(xyz1, ijk1(:n), end_in_mesh)
! If this is the first iteration of the filter loop, check if the track
! intersects any part of the mesh.
if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then
if (.not. m % intersects(xyz0, xyz1)) return
end if
! ! If this is the first iteration of the filter loop, check if the track
! ! intersects any part of the mesh.
! if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then
! if (.not. m % intersects(xyz0, xyz1)) return
! end if
! ========================================================================
! Figure out which mesh cell to tally.
! ! ========================================================================
! ! Figure out which mesh cell to tally.
! Copy the un-modified coordinates the particle direction.
xyz0 = p % last_xyz
xyz1 = p % coord(1) % xyz
uvw = p % coord(1) % uvw
! ! Copy the un-modified coordinates the particle direction.
! xyz0 = p % last_xyz
! xyz1 = p % coord(1) % xyz
! uvw = p % coord(1) % uvw
! Compute the length of the entire track.
total_distance = sqrt(sum((xyz1 - xyz0)**2))
! ! Compute the length of the entire track.
! total_distance = sqrt(sum((xyz1 - xyz0)**2))
! We are looking for the first valid mesh bin. Check to see if the
! particle starts inside the mesh.
if (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension)) then
! The particle does not start in the mesh. Note that we nudged the
! start and end coordinates by a TINY_BIT each so we will have
! difficulty resolving tracks that are less than 2*TINY_BIT in length.
! If the track is that short, it is also insignificant so we can
! safely ignore it in the tallies.
if (total_distance < 2*TINY_BIT) return
! ! We are looking for the first valid mesh bin. Check to see if the
! ! particle starts inside the mesh.
! if (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension)) then
! ! The particle does not start in the mesh. Note that we nudged the
! ! start and end coordinates by a TINY_BIT each so we will have
! ! difficulty resolving tracks that are less than 2*TINY_BIT in length.
! ! If the track is that short, it is also insignificant so we can
! ! safely ignore it in the tallies.
! if (total_distance < 2*TINY_BIT) return
! The particle does not start in the mesh so keep iterating the ijk0
! indices to cross the nearest mesh surface until we've found a valid
! bin. MAX_SEARCH_ITER prevents an infinite loop.
search_iter = 0
do while (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension))
if (search_iter == MAX_SEARCH_ITER) then
call warning("Failed to find a mesh intersection on a tally mesh &
&filter.")
return
end if
! ! The particle does not start in the mesh so keep iterating the ijk0
! ! indices to cross the nearest mesh surface until we've found a valid
! ! bin. MAX_SEARCH_ITER prevents an infinite loop.
! search_iter = 0
! do while (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension))
! if (search_iter == MAX_SEARCH_ITER) then
! call warning("Failed to find a mesh intersection on a tally mesh &
! &filter.")
! return
! end if
do j = 1, n
if (abs(uvw(j)) < FP_PRECISION) then
d(j) = INFINITY
else if (uvw(j) > 0) then
xyz_cross = m % lower_left(j) + ijk0(j) * m % width(j)
d(j) = (xyz_cross - xyz0(j)) / uvw(j)
else
xyz_cross = m % lower_left(j) + (ijk0(j) - 1) * m % width(j)
d(j) = (xyz_cross - xyz0(j)) / uvw(j)
end if
end do
j = minloc(d(:n), 1)
if (uvw(j) > ZERO) then
ijk0(j) = ijk0(j) + 1
else
ijk0(j) = ijk0(j) - 1
end if
! do j = 1, n
! if (abs(uvw(j)) < FP_PRECISION) then
! d(j) = INFINITY
! else if (uvw(j) > 0) then
! xyz_cross = m % lower_left(j) + ijk0(j) * m % width(j)
! d(j) = (xyz_cross - xyz0(j)) / uvw(j)
! else
! xyz_cross = m % lower_left(j) + (ijk0(j) - 1) * m % width(j)
! d(j) = (xyz_cross - xyz0(j)) / uvw(j)
! end if
! end do
! j = minloc(d(:n), 1)
! if (uvw(j) > ZERO) then
! ijk0(j) = ijk0(j) + 1
! else
! ijk0(j) = ijk0(j) - 1
! end if
search_iter = search_iter + 1
end do
distance = d(j)
xyz0 = xyz0 + distance * uvw
end if
! search_iter = search_iter + 1
! end do
! distance = d(j)
! xyz0 = xyz0 + distance * uvw
! end if
do
! ========================================================================
! Compute the length of the track segment in the appropiate mesh cell and
! return.
! do
! ! ========================================================================
! ! Compute the length of the track segment in the appropiate mesh cell and
! ! return.
if (all(ijk0(:n) == ijk1(:n))) then
! The track ends in this cell. Use the particle end location rather
! than the mesh surface.
distance = sqrt(sum((xyz1 - xyz0)**2))
else
! The track exits this cell. Determine the distance to the closest mesh
! surface.
do j = 1, n
if (abs(uvw(j)) < FP_PRECISION) then
d(j) = INFINITY
else if (uvw(j) > 0) then
xyz_cross = m % lower_left(j) + ijk0(j) * m % width(j)
d(j) = (xyz_cross - xyz0(j)) / uvw(j)
else
xyz_cross = m % lower_left(j) + (ijk0(j) - 1) * m % width(j)
d(j) = (xyz_cross - xyz0(j)) / uvw(j)
end if
end do
j = minloc(d(:n), 1)
distance = d(j)
end if
! if (all(ijk0(:n) == ijk1(:n))) then
! ! The track ends in this cell. Use the particle end location rather
! ! than the mesh surface.
! distance = sqrt(sum((xyz1 - xyz0)**2))
! else
! ! The track exits this cell. Determine the distance to the closest mesh
! ! surface.
! do j = 1, n
! if (abs(uvw(j)) < FP_PRECISION) then
! d(j) = INFINITY
! else if (uvw(j) > 0) then
! xyz_cross = m % lower_left(j) + ijk0(j) * m % width(j)
! d(j) = (xyz_cross - xyz0(j)) / uvw(j)
! else
! xyz_cross = m % lower_left(j) + (ijk0(j) - 1) * m % width(j)
! d(j) = (xyz_cross - xyz0(j)) / uvw(j)
! end if
! end do
! j = minloc(d(:n), 1)
! distance = d(j)
! end if
! Assign the next tally bin and the score.
bin = m % get_bin_from_indices(ijk0(:n))
call match % bins % push_back(bin)
call match % weights % push_back(distance / total_distance)
! ! Assign the next tally bin and the score.
! bin = m % get_bin_from_indices(ijk0(:n))
! call match % bins % push_back(bin)
! call match % weights % push_back(distance / total_distance)
! Find the next mesh cell that the particle enters.
! ! Find the next mesh cell that the particle enters.
! If the particle track ends in that bin, then we are done.
if (all(ijk0(:n) == ijk1(:n))) exit
! ! If the particle track ends in that bin, then we are done.
! if (all(ijk0(:n) == ijk1(:n))) exit
! Translate the starting coordintes by the distance to that face. This
! should be the xyz that we computed the distance to in the last
! iteration of the filter loop.
xyz0 = xyz0 + distance * uvw
! ! Translate the starting coordintes by the distance to that face. This
! ! should be the xyz that we computed the distance to in the last
! ! iteration of the filter loop.
! xyz0 = xyz0 + distance * uvw
! Increment the indices into the next mesh cell.
if (uvw(j) > ZERO) then
ijk0(j) = ijk0(j) + 1
else
ijk0(j) = ijk0(j) - 1
end if
! ! Increment the indices into the next mesh cell.
! if (uvw(j) > ZERO) then
! ijk0(j) = ijk0(j) + 1
! else
! ijk0(j) = ijk0(j) - 1
! end if
! If the next indices are invalid, then the track has left the mesh and
! we are done.
if (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension)) exit
end do
! ! If the next indices are invalid, then the track has left the mesh and
! ! we are done.
! if (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension)) exit
! end do
end subroutine get_all_bins_mesh
@ -248,9 +247,12 @@ contains
class(MeshFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
type(RegularMesh) :: m
m = meshes(this % mesh)
call write_dataset(filter_group, "type", "mesh")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "bins", meshes(this % mesh) % id)
call write_dataset(filter_group, "bins", m % id())
end subroutine to_statepoint_mesh
function text_label_mesh(this, bin) result(label)
@ -259,20 +261,20 @@ contains
character(MAX_LINE_LEN) :: label
integer, allocatable :: ijk(:)
type(RegularMesh) :: m
associate (m => meshes(this % mesh))
allocate(ijk(m % n_dimension))
call m % get_indices_from_bin(bin, ijk)
if (m % n_dimension == 1) then
label = "Mesh Index (" // trim(to_str(ijk(1))) // ")"
elseif (m % n_dimension == 2) then
label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ")"
elseif (m % n_dimension == 3) then
label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")"
end if
end associate
m = meshes(this % mesh)
allocate(ijk(m % n_dimension()))
call m % get_indices_from_bin(bin, ijk)
if (m % n_dimension() == 1) then
label = "Mesh Index (" // trim(to_str(ijk(1))) // ")"
elseif (m % n_dimension() == 2) then
label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ")"
elseif (m % n_dimension() == 3) then
label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")"
end if
end function text_label_mesh
!===============================================================================
@ -304,18 +306,25 @@ contains
integer(C_INT32_T), value, intent(in) :: index_mesh
integer(C_INT) :: err
type(RegularMesh) :: m
integer :: i
err = verify_filter(index)
if (err == 0) then
select type (f => filters(index) % obj)
type is (MeshFilter)
if (index_mesh >= 1 .and. index_mesh <= n_meshes) then
if (index_mesh >= 0 .and. index_mesh < n_meshes()) then
f % mesh = index_mesh
f % n_bins = product(meshes(index_mesh) % dimension)
f % n_bins = 1
m = meshes(index_mesh)
do i = 1, m % n_dimension()
f % n_bins = f % n_bins * m % dimension(i)
end do
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in 'meshes' array is out of bounds.")
end if
class default
class default
err = E_INVALID_TYPE
call set_errmsg("Tried to set mesh on a non-mesh filter.")
end select

309
src/tallies/tally_filter_meshsurface.F90
View file

@ -42,6 +42,7 @@ contains
integer :: id
integer :: n
integer :: n_dim
integer(C_INT) :: err
type(RegularMesh) :: m
n = node_word_count(node, "bins")
@ -90,133 +91,130 @@ contains
real(8) :: distance ! actual distance traveled
logical :: start_in_mesh ! particle's starting xyz in mesh?
logical :: end_in_mesh ! particle's ending xyz in mesh?
type(RegularMesh) :: m
! Copy starting and ending location of particle
xyz0 = p % last_xyz_current
xyz1 = p % coord(1) % xyz
! ! Copy starting and ending location of particle
! xyz0 = p % last_xyz_current
! xyz1 = p % coord(1) % xyz
associate (m => meshes(this % mesh))
n_dim = m % n_dimension
! m = meshes(this % mesh)
! n_dim = m % n_dimension()
! Determine indices for starting and ending location
call m % get_indices(xyz0, ijk0, start_in_mesh)
call m % get_indices(xyz1, ijk1, end_in_mesh)
! ! Determine indices for starting and ending location
! call m % get_indices(xyz0, ijk0, start_in_mesh)
! call m % get_indices(xyz1, ijk1, end_in_mesh)
! Check to see if start or end is in mesh -- if not, check if track still
! intersects with mesh
if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then
if (.not. m % intersects(xyz0, xyz1)) return
end if
! ! Check to see if start or end is in mesh -- if not, check if track still
! ! intersects with mesh
! if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then
! !if (.not. m % intersects(xyz0, xyz1)) return
! end if
! Calculate number of surface crossings
n_cross = sum(abs(ijk1(:n_dim) - ijk0(:n_dim)))
if (n_cross == 0) return
! ! Calculate number of surface crossings
! n_cross = sum(abs(ijk1(:n_dim) - ijk0(:n_dim)))
! if (n_cross == 0) return
! Copy particle's direction
uvw = p % coord(1) % uvw
! ! Copy particle's direction
! uvw = p % coord(1) % uvw
! Bounding coordinates
do d1 = 1, n_dim
if (uvw(d1) > 0) then
xyz_cross(d1) = m % lower_left(d1) + ijk0(d1) * m % width(d1)
else
xyz_cross(d1) = m % lower_left(d1) + (ijk0(d1) - 1) * m % width(d1)
end if
end do
! ! Bounding coordinates
! do d1 = 1, n_dim
! if (uvw(d1) > 0) then
! xyz_cross(d1) = m % lower_left(d1) + ijk0(d1) * m % width(d1)
! else
! xyz_cross(d1) = m % lower_left(d1) + (ijk0(d1) - 1) * m % width(d1)
! end if
! end do
do j = 1, n_cross
! Set the distances to infinity
d = INFINITY
! do j = 1, n_cross
! ! Set the distances to infinity
! d = INFINITY
! Calculate distance to each bounding surface. We need to treat
! special case where the cosine of the angle is zero since this would
! result in a divide-by-zero.
do d1 = 1, n_dim
if (uvw(d1) == 0) then
d(d1) = INFINITY
else
d(d1) = (xyz_cross(d1) - xyz0(d1))/uvw(d1)
end if
end do
! ! Calculate distance to each bounding surface. We need to treat
! ! special case where the cosine of the angle is zero since this would
! ! result in a divide-by-zero.
! do d1 = 1, n_dim
! if (uvw(d1) == 0) then
! d(d1) = INFINITY
! else
! d(d1) = (xyz_cross(d1) - xyz0(d1))/uvw(d1)
! end if
! end do
! Determine the closest bounding surface of the mesh cell by
! calculating the minimum distance. Then use the minimum distance and
! direction of the particle to determine which surface was crossed.
distance = minval(d)
! ! Determine the closest bounding surface of the mesh cell by
! ! calculating the minimum distance. Then use the minimum distance and
! ! direction of the particle to determine which surface was crossed.
! distance = minval(d)
! Loop over the dimensions
do d1 = 1, n_dim
! ! Loop over the dimensions
! do d1 = 1, n_dim
! Check whether distance is the shortest distance
if (distance == d(d1)) then
! ! Check whether distance is the shortest distance
! if (distance == d(d1)) then
! Check whether particle is moving in positive d1 direction
if (uvw(d1) > 0) then
! ! Check whether particle is moving in positive d1 direction
! if (uvw(d1) > 0) then
! Outward current on d1 max surface
if (all(ijk0(:n_dim) >= 1) .and. &
all(ijk0(:n_dim) <= m % dimension)) then
i_surf = d1 * 4 - 1
i_mesh = m % get_bin_from_indices(ijk0)
i_bin = 4*n_dim*(i_mesh - 1) + i_surf
! ! Outward current on d1 max surface
! if (all(ijk0(:n_dim) >= 1) .and. &
! all(ijk0(:n_dim) <= m % dimension)) then
! i_surf = d1 * 4 - 1
! i_mesh = m % get_bin_from_indices(ijk0)
! i_bin = 4*n_dim*(i_mesh - 1) + i_surf
call match % bins % push_back(i_bin)
call match % weights % push_back(ONE)
end if
! call match % bins % push_back(i_bin)
! call match % weights % push_back(ONE)
! end if
! Advance position
ijk0(d1) = ijk0(d1) + 1
xyz_cross(d1) = xyz_cross(d1) + m % width(d1)
! ! Advance position
! ijk0(d1) = ijk0(d1) + 1
! xyz_cross(d1) = xyz_cross(d1) + m % width(d1)
! If the particle crossed the surface, tally the inward current on
! d1 min surface
if (all(ijk0(:n_dim) >= 1) .and. &
all(ijk0(:n_dim) <= m % dimension)) then
i_surf = d1 * 4 - 2
i_mesh = m % get_bin_from_indices(ijk0)
i_bin = 4*n_dim*(i_mesh - 1) + i_surf
! ! If the particle crossed the surface, tally the inward current on
! ! d1 min surface
! if (all(ijk0(:n_dim) >= 1)) then ! .and. all(ijk0(:n_dim) <= m % dimension)) then
! i_surf = d1 * 4 - 2
! i_mesh = m % get_bin_from_indices(ijk0)
! i_bin = 4*n_dim*(i_mesh - 1) + i_surf
call match % bins % push_back(i_bin)
call match % weights % push_back(ONE)
end if
! call match % bins % push_back(i_bin)
! call match % weights % push_back(ONE)
! end if
else
! The particle is moving in the negative d1 direction
! else
! ! The particle is moving in the negative d1 direction
! Outward current on d1 min surface
if (all(ijk0(:n_dim) >= 1) .and. &
all(ijk0(:n_dim) <= m % dimension)) then
i_surf = d1 * 4 - 3
i_mesh = m % get_bin_from_indices(ijk0)
i_bin = 4*n_dim*(i_mesh - 1) + i_surf
! ! Outward current on d1 min surface
! if (all(ijk0(:n_dim) >= 1)) then ! .and. all(ijk0(:n_dim) <= m % dimension)) then
! i_surf = d1 * 4 - 3
! i_mesh = m % get_bin_from_indices(ijk0)
! i_bin = 4*n_dim*(i_mesh - 1) + i_surf
call match % bins % push_back(i_bin)
call match % weights % push_back(ONE)
end if
! call match % bins % push_back(i_bin)
! call match % weights % push_back(ONE)
! end if
! Advance position
ijk0(d1) = ijk0(d1) - 1
xyz_cross(d1) = xyz_cross(d1) - m % width(d1)
! ! Advance position
! ijk0(d1) = ijk0(d1) - 1
! xyz_cross(d1) = xyz_cross(d1) - m % width(d1)
! If the particle crossed the surface, tally the inward current on
! d1 max surface
if (all(ijk0(:n_dim) >= 1) .and. &
all(ijk0(:n_dim) <= m % dimension)) then
i_surf = d1 * 4
i_mesh = m % get_bin_from_indices(ijk0)
i_bin = 4*n_dim*(i_mesh - 1) + i_surf
! ! If the particle crossed the surface, tally the inward current on
! ! d1 max surface
! if (all(ijk0(:n_dim) >= 1)) then ! .and. all(ijk0(:n_dim) <= m % dimension)) then
! i_surf = d1 * 4
! i_mesh = m % get_bin_from_indices(ijk0)
! i_bin = 4*n_dim*(i_mesh - 1) + i_surf
call match % bins % push_back(i_bin)
call match % weights % push_back(ONE)
end if
end if
end if
end do
! call match % bins % push_back(i_bin)
! call match % weights % push_back(ONE)
! end if
! end if
! end if
! end do
! Calculate new coordinates
xyz0 = xyz0 + distance * uvw
end do
end associate
! ! Calculate new coordinates
! xyz0 = xyz0 + distance * uvw
! end do
end subroutine get_all_bins
@ -224,9 +222,12 @@ contains
class(MeshSurfaceFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
type(RegularMesh) :: m
m = meshes(this % mesh)
call write_dataset(filter_group, "type", "meshsurface")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "bins", meshes(this % mesh) % id)
call write_dataset(filter_group, "bins", m % id())
end subroutine to_statepoint
function text_label(this, bin) result(label)
@ -238,55 +239,55 @@ contains
integer :: i_surf
integer :: n_dim
integer, allocatable :: ijk(:)
type(RegularMesh) :: m
associate (m => meshes(this % mesh))
n_dim = m % n_dimension
allocate(ijk(n_dim))
m = meshes(this % mesh)
n_dim = m % n_dimension()
allocate(ijk(n_dim))
! Get flattend mesh index and surface index
i_mesh = (bin - 1) / (4*n_dim) + 1
i_surf = mod(bin - 1, 4*n_dim) + 1
! Get flattend mesh index and surface index
i_mesh = (bin - 1) / (4*n_dim) + 1
i_surf = mod(bin - 1, 4*n_dim) + 1
! Get mesh index part of label
call m % get_indices_from_bin(i_mesh, ijk)
if (m % n_dimension == 1) then
label = "Mesh Index (" // trim(to_str(ijk(1))) // ")"
elseif (m % n_dimension == 2) then
label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ")"
elseif (m % n_dimension == 3) then
label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")"
end if
! Get mesh index part of label
call m % get_indices_from_bin(i_mesh, ijk)
if (m % n_dimension() == 1) then
label = "Mesh Index (" // trim(to_str(ijk(1))) // ")"
elseif (m % n_dimension() == 2) then
label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ")"
elseif (m % n_dimension() == 3) then
label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")"
end if
! Get surface part of label
select case (i_surf)
case (OUT_LEFT)
label = trim(label) // " Outgoing, x-min"
case (IN_LEFT)
label = trim(label) // " Incoming, x-min"
case (OUT_RIGHT)
label = trim(label) // " Outgoing, x-max"
case (IN_RIGHT)
label = trim(label) // " Incoming, x-max"
case (OUT_BACK)
label = trim(label) // " Outgoing, y-min"
case (IN_BACK)
label = trim(label) // " Incoming, y-min"
case (OUT_FRONT)
label = trim(label) // " Outgoing, y-max"
case (IN_FRONT)
label = trim(label) // " Incoming, y-max"
case (OUT_BOTTOM)
label = trim(label) // " Outgoing, z-min"
case (IN_BOTTOM)
label = trim(label) // " Incoming, z-min"
case (OUT_TOP)
label = trim(label) // " Outgoing, z-max"
case (IN_TOP)
label = trim(label) // " Incoming, z-max"
end select
end associate
! Get surface part of label
select case (i_surf)
case (OUT_LEFT)
label = trim(label) // " Outgoing, x-min"
case (IN_LEFT)
label = trim(label) // " Incoming, x-min"
case (OUT_RIGHT)
label = trim(label) // " Outgoing, x-max"
case (IN_RIGHT)
label = trim(label) // " Incoming, x-max"
case (OUT_BACK)
label = trim(label) // " Outgoing, y-min"
case (IN_BACK)
label = trim(label) // " Incoming, y-min"
case (OUT_FRONT)
label = trim(label) // " Outgoing, y-max"
case (IN_FRONT)
label = trim(label) // " Incoming, y-max"
case (OUT_BOTTOM)
label = trim(label) // " Outgoing, z-min"
case (IN_BOTTOM)
label = trim(label) // " Incoming, z-min"
case (OUT_TOP)
label = trim(label) // " Outgoing, z-max"
case (IN_TOP)
label = trim(label) // " Incoming, z-max"
end select
end function text_label
!===============================================================================
@ -318,16 +319,22 @@ contains
integer(C_INT32_T), value, intent(in) :: index_mesh
integer(C_INT) :: err
integer :: i
integer :: n_dim
type(RegularMesh) :: m
err = verify_filter(index)
if (err == 0) then
select type (f => filters(index) % obj)
type is (MeshSurfaceFilter)
if (index_mesh >= 1 .and. index_mesh <= n_meshes) then
if (index_mesh >= 0 .and. index_mesh < n_meshes()) then
f % mesh = index_mesh
n_dim = meshes(index_mesh) % n_dimension
f % n_bins = 4*n_dim*product(meshes(index_mesh) % dimension)
m = meshes(index_mesh)
n_dim = m % n_dimension()
f % n_bins = 4*n_dim
do i = 1, n_dim
f % n_bins = f % n_bins * m % dimension(i)
end do
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in 'meshes' array is out of bounds.")

11
src/tallies/trigger.F90
View file

@ -257,14 +257,14 @@ contains
logical :: print_ebin ! should incoming energy bin be displayed?
real(8) :: rel_err = ZERO ! temporary relative error of result
real(8) :: std_dev = ZERO ! temporary standard deviration of result
type(RegularMesh), pointer :: m ! surface current mesh
type(RegularMesh) :: m ! surface current mesh
! Get pointer to mesh
i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
i_filter_surf = t % filter(t % find_filter(FILTER_SURFACE))
select type(filt => filters(i_filter_mesh) % obj)
type is (MeshFilter)
m => meshes(filt % mesh)
m = meshes(filt % mesh)
end select
! initialize bins array
@ -285,8 +285,11 @@ contains
end if
! Get the dimensions and number of cells in the mesh
n_dim = m % n_dimension
n_cells = product(m % dimension)
n_dim = m % n_dimension()
n_cells = 1
do j = 1, n_dim
n_cells = n_cells * m % dimension(j)
end do
! Loop over all the mesh cells
do i = 1, n_cells