diff --git a/CMakeLists.txt b/CMakeLists.txt index a13f30b575..ec77fa0e33 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -431,10 +431,12 @@ set(LIBOPENMC_FORTRAN_SRC src/tallies/trigger_header.F90 ) set(LIBOPENMC_CXX_SRC - src/constants.h + src/cell.cpp src/initialize.cpp src/finalize.cpp + src/geometry_aux.cpp src/hdf5_interface.cpp + src/lattice.cpp src/math_functions.cpp src/message_passing.cpp src/mgxs.cpp diff --git a/docs/source/devguide/styleguide.rst b/docs/source/devguide/styleguide.rst index 89d5673943..ccf6bf5217 100644 --- a/docs/source/devguide/styleguide.rst +++ b/docs/source/devguide/styleguide.rst @@ -219,8 +219,8 @@ Curly braces For a function definition, the opening and closing braces should each be on their own lines. This helps distinguish function code from the argument list. -If the entire function fits on one line, then the braces can be on the same -line. e.g.: +If the entire function fits on one or two lines, then the braces can be on the +same line. e.g.: .. code-block:: C++ @@ -238,6 +238,9 @@ line. e.g.: int return_one() {return 1;} + int return_one() + {return 1;} + For a conditional, the opening brace should be on the same line as the end of the conditional statement. If there is a following ``else if`` or ``else`` statement, the closing brace should be on the same line as that following diff --git a/examples/python/pincell_depletion/restart_depletion.py b/examples/python/pincell_depletion/restart_depletion.py new file mode 100644 index 0000000000..f9324a8cb8 --- /dev/null +++ b/examples/python/pincell_depletion/restart_depletion.py @@ -0,0 +1,103 @@ +import openmc +import openmc.deplete +import numpy as np +import matplotlib.pyplot as plt + +############################################################################### +# Simulation Input File Parameters +############################################################################### + +# OpenMC simulation parameters +batches = 100 +inactive = 10 +particles = 1000 + +# Depletion simulation parameters +time_step = 1*24*60*60 # s +final_time = 5*24*60*60 # s +time_steps = np.full(final_time // time_step, time_step) + +chain_file = './chain_simple.xml' +power = 174 # W/cm, for 2D simulations only (use W for 3D) + +############################################################################### +# Load previous simulation results +############################################################################### + +# Load geometry from statepoint +statepoint = 'statepoint.100.h5' +with openmc.StatePoint(statepoint) as sp: + geometry = sp.summary.geometry + +# Load previous depletion results +previous_results = openmc.deplete.ResultsList("depletion_results.h5") + +############################################################################### +# Transport calculation settings +############################################################################### + +# Instantiate a Settings object, set all runtime parameters +settings_file = openmc.Settings() +settings_file.batches = batches +settings_file.inactive = inactive +settings_file.particles = particles + +# Create an initial uniform spatial source distribution over fissionable zones +bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] +uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True) +settings_file.source = openmc.source.Source(space=uniform_dist) + +entropy_mesh = openmc.Mesh() +entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50] +entropy_mesh.upper_right = [0.39218, 0.39218, 1.e50] +entropy_mesh.dimension = [10, 10, 1] +settings_file.entropy_mesh = entropy_mesh + +############################################################################### +# Initialize and run depletion calculation +############################################################################### + +op = openmc.deplete.Operator(geometry, settings_file, chain_file, + previous_results) + +# Perform simulation using the predictor algorithm +openmc.deplete.integrator.predictor(op, time_steps, power) + +############################################################################### +# Read depletion calculation results +############################################################################### + +# Open results file +results = openmc.deplete.ResultsList("depletion_results.h5") + +# Obtain K_eff as a function of time +time, keff = results.get_eigenvalue() + +# Obtain U235 concentration as a function of time +time, n_U235 = results.get_atoms('1', 'U235') + +# Obtain Xe135 absorption as a function of time +time, Xe_gam = results.get_reaction_rate('1', 'Xe135', '(n,gamma)') + +############################################################################### +# Generate plots +############################################################################### + +plt.figure() +plt.plot(time/(24*60*60), keff, label="K-effective") +plt.xlabel("Time (days)") +plt.ylabel("Keff") +plt.show() + +plt.figure() +plt.plot(time/(24*60*60), n_U235, label="U 235") +plt.xlabel("Time (days)") +plt.ylabel("n U5 (-)") +plt.show() + +plt.figure() +plt.plot(time/(24*60*60), Xe_gam, label="Xe135 absorption") +plt.xlabel("Time (days)") +plt.ylabel("RR (-)") +plt.show() +plt.close('all') diff --git a/examples/python/pincell_depletion/run_depletion.py b/examples/python/pincell_depletion/run_depletion.py index 8c12211bca..a2a5e0e7f2 100644 --- a/examples/python/pincell_depletion/run_depletion.py +++ b/examples/python/pincell_depletion/run_depletion.py @@ -1,6 +1,7 @@ import openmc import openmc.deplete import numpy as np +import matplotlib.pyplot as plt ############################################################################### # Simulation Input File Parameters @@ -13,7 +14,7 @@ particles = 1000 # Depletion simulation parameters time_step = 1*24*60*60 # s -final_time = 15*24*60*60 # s +final_time = 5*24*60*60 # s time_steps = np.full(final_time // time_step, time_step) chain_file = './chain_simple.xml' @@ -49,7 +50,7 @@ borated_water.add_element('O', 2.4e-2) borated_water.add_s_alpha_beta('c_H_in_H2O') ############################################################################### -# Exporting to OpenMC geometry.xml file +# Create geometry ############################################################################### # Instantiate ZCylinder surfaces @@ -94,7 +95,7 @@ root.add_cells([fuel, gap, clad, water]) geometry = openmc.Geometry(root) ############################################################################### -# Exporting to OpenMC materials.xml file +# Set volumes of depletable materials ############################################################################### # Compute cell areas @@ -105,7 +106,7 @@ area[fuel] = np.pi * fuel_or.coefficients['R'] ** 2 uo2.volume = area[fuel] ############################################################################### -# Exporting to OpenMC settings.xml file +# Transport calculation settings ############################################################################### # Instantiate a Settings object, set all runtime parameters, and export to XML @@ -143,6 +144,32 @@ results = openmc.deplete.ResultsList("depletion_results.h5") # Obtain K_eff as a function of time time, keff = results.get_eigenvalue() - + # Obtain U235 concentration as a function of time time, n_U235 = results.get_atoms('1', 'U235') + +# Obtain Xe135 absorption as a function of time +time, Xe_gam = results.get_reaction_rate('1', 'Xe135', '(n,gamma)') + +############################################################################### +# Generate plots +############################################################################### + +plt.figure() +plt.plot(time/(24*60*60), keff, label="K-effective") +plt.xlabel("Time (days)") +plt.ylabel("Keff") +plt.show() + +plt.figure() +plt.plot(time/(24*60*60), n_U235, label="U 235") +plt.xlabel("Time (days)") +plt.ylabel("n U5 (-)") +plt.show() + +plt.figure() +plt.plot(time/(24*60*60), Xe_gam, label="Xe135 absorption") +plt.xlabel("Time (days)") +plt.ylabel("RR (-)") +plt.show() +plt.close('all') diff --git a/include/openmc.h b/include/openmc.h index 38b7f21670..2c0f10eae0 100644 --- a/include/openmc.h +++ b/include/openmc.h @@ -17,6 +17,7 @@ extern "C" { }; int openmc_calculate_volumes(); + int openmc_cell_filter_get_bins(int32_t index, int32_t** cells, int32_t* n); int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n); int openmc_cell_get_id(int32_t index, int32_t* id); int openmc_cell_set_fill(int32_t index, int type, int32_t n, const int32_t* indices); @@ -31,7 +32,7 @@ extern "C" { int openmc_extend_sources(int32_t n, int32_t* index_start, int32_t* index_end); int openmc_extend_tallies(int32_t n, int32_t* index_start, int32_t* index_end); int openmc_filter_get_id(int32_t index, int32_t* id); - int openmc_filter_get_type(int32_t index, const char** type); + int openmc_filter_get_type(int32_t index, char* type); int openmc_filter_set_id(int32_t index, int32_t id); int openmc_filter_set_type(int32_t index, const char* type); int openmc_finalize(); diff --git a/openmc/capi/filter.py b/openmc/capi/filter.py index 5c21c3ef4b..391e63f073 100644 --- a/openmc/capi/filter.py +++ b/openmc/capi/filter.py @@ -17,11 +17,16 @@ from .mesh import Mesh __all__ = ['Filter', 'AzimuthalFilter', 'CellFilter', 'CellbornFilter', 'CellfromFilter', 'DistribcellFilter', 'DelayedGroupFilter', 'EnergyFilter', 'EnergyoutFilter', - 'EnergyFunctionFilter', 'MaterialFilter', 'MeshFilter', - 'MeshSurfaceFilter', 'MuFilter', 'PolarFilter', 'SurfaceFilter', - 'UniverseFilter', 'filters'] + 'EnergyFunctionFilter', 'LegendreFilter', 'MaterialFilter', 'MeshFilter', + 'MeshSurfaceFilter', 'MuFilter', 'PolarFilter', 'SphericalHarmonicsFilter', + 'SpatialLegendreFilter', 'SurfaceFilter', + 'UniverseFilter', 'ZernikeFilter', 'filters'] # Tally functions +_dll.openmc_cell_filter_get_bins.argtypes = [ + c_int32, POINTER(POINTER(c_int32)), POINTER(c_int32)] +_dll.openmc_cell_filter_get_bins.restype = c_int +_dll.openmc_cell_filter_get_bins.errcheck = _error_handler _dll.openmc_energy_filter_get_bins.argtypes = [ c_int32, POINTER(POINTER(c_double)), POINTER(c_int32)] _dll.openmc_energy_filter_get_bins.restype = c_int @@ -47,6 +52,12 @@ _dll.openmc_filter_set_type.errcheck = _error_handler _dll.openmc_get_filter_index.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_get_filter_index.restype = c_int _dll.openmc_get_filter_index.errcheck = _error_handler +_dll.openmc_legendre_filter_get_order.argtypes = [c_int32, POINTER(c_int)] +_dll.openmc_legendre_filter_get_order.restype = c_int +_dll.openmc_legendre_filter_get_order.errcheck = _error_handler +_dll.openmc_legendre_filter_set_order.argtypes = [c_int32, c_int] +_dll.openmc_legendre_filter_set_order.restype = c_int +_dll.openmc_legendre_filter_set_order.errcheck = _error_handler _dll.openmc_material_filter_get_bins.argtypes = [ c_int32, POINTER(POINTER(c_int32)), POINTER(c_int32)] _dll.openmc_material_filter_get_bins.restype = c_int @@ -66,7 +77,24 @@ _dll.openmc_meshsurface_filter_get_mesh.errcheck = _error_handler _dll.openmc_meshsurface_filter_set_mesh.argtypes = [c_int32, c_int32] _dll.openmc_meshsurface_filter_set_mesh.restype = c_int _dll.openmc_meshsurface_filter_set_mesh.errcheck = _error_handler - +_dll.openmc_spatial_legendre_filter_get_order.argtypes = [c_int32, POINTER(c_int)] +_dll.openmc_spatial_legendre_filter_get_order.restype = c_int +_dll.openmc_spatial_legendre_filter_get_order.errcheck = _error_handler +_dll.openmc_spatial_legendre_filter_set_order.argtypes = [c_int32, c_int] +_dll.openmc_spatial_legendre_filter_set_order.restype = c_int +_dll.openmc_spatial_legendre_filter_set_order.errcheck = _error_handler +_dll.openmc_sphharm_filter_get_order.argtypes = [c_int32, POINTER(c_int)] +_dll.openmc_sphharm_filter_get_order.restype = c_int +_dll.openmc_sphharm_filter_get_order.errcheck = _error_handler +_dll.openmc_sphharm_filter_set_order.argtypes = [c_int32, c_int] +_dll.openmc_sphharm_filter_set_order.restype = c_int +_dll.openmc_sphharm_filter_set_order.errcheck = _error_handler +_dll.openmc_zernike_filter_get_order.argtypes = [c_int32, POINTER(c_int)] +_dll.openmc_zernike_filter_get_order.restype = c_int +_dll.openmc_zernike_filter_get_order.errcheck = _error_handler +_dll.openmc_zernike_filter_set_order.argtypes = [c_int32, c_int] +_dll.openmc_zernike_filter_set_order.restype = c_int +_dll.openmc_zernike_filter_set_order.errcheck = _error_handler class Filter(_FortranObjectWithID): __instances = WeakValueDictionary() @@ -150,6 +178,13 @@ class AzimuthalFilter(Filter): class CellFilter(Filter): filter_type = 'cell' + @property + def bins(self): + cells = POINTER(c_int32)() + n = c_int32() + _dll.openmc_cell_filter_get_bins(self._index, cells, n) + return as_array(cells, (n.value,)) + class CellbornFilter(Filter): filter_type = 'cellborn' @@ -171,6 +206,25 @@ class EnergyFunctionFilter(Filter): filter_type = 'energyfunction' +class LegendreFilter(Filter): + filter_type = 'legendre' + + def __init__(self, order=None, uid=None, new=True, index=None): + super().__init__(uid, new, index) + if order is not None: + self.order = order + + @property + def order(self): + temp_order = c_int() + _dll.openmc_legendre_filter_get_order(self._index, temp_order) + return temp_order.value + + @order.setter + def order(self, order): + _dll.openmc_legendre_filter_set_order(self._index, order) + + class MaterialFilter(Filter): filter_type = 'material' @@ -241,6 +295,44 @@ class PolarFilter(Filter): filter_type = 'polar' +class SphericalHarmonicsFilter(Filter): + filter_type = 'sphericalharmonics' + + def __init__(self, order=None, uid=None, new=True, index=None): + super().__init__(uid, new, index) + if order is not None: + self.order = order + + @property + def order(self): + temp_order = c_int() + _dll.openmc_sphharm_filter_get_order(self._index, temp_order) + return temp_order.value + + @order.setter + def order(self, order): + _dll.openmc_sphharm_filter_set_order(self._index, order) + + +class SpatialLegendreFilter(Filter): + filter_type = 'spatiallegendre' + + def __init__(self, order=None, uid=None, new=True, index=None): + super().__init__(uid, new, index) + if order is not None: + self.order = order + + @property + def order(self): + temp_order = c_int() + _dll.openmc_spatial_legendre_filter_get_order(self._index, temp_order) + return temp_order.value + + @order.setter + def order(self, order): + _dll.openmc_spatial_legendre_filter_set_order(self._index, order) + + class SurfaceFilter(Filter): filter_type = 'surface' @@ -249,6 +341,25 @@ class UniverseFilter(Filter): filter_type = 'universe' +class ZernikeFilter(Filter): + filter_type = 'zernike' + + def __init__(self, order=None, uid=None, new=True, index=None): + super().__init__(uid, new, index) + if order is not None: + self.order = order + + @property + def order(self): + temp_order = c_int() + _dll.openmc_zernike_filter_get_order(self._index, temp_order) + return temp_order.value + + @order.setter + def order(self, order): + _dll.openmc_zernike_filter_set_order(self._index, order) + + _FILTER_TYPE_MAP = { 'azimuthal': AzimuthalFilter, 'cell': CellFilter, @@ -259,13 +370,17 @@ _FILTER_TYPE_MAP = { 'energy': EnergyFilter, 'energyout': EnergyoutFilter, 'energyfunction': EnergyFunctionFilter, + 'legendre': LegendreFilter, 'material': MaterialFilter, 'mesh': MeshFilter, 'meshsurface': MeshSurfaceFilter, 'mu': MuFilter, 'polar': PolarFilter, + 'sphericalharmonics': SphericalHarmonicsFilter, + 'spatiallegendre': SpatialLegendreFilter, 'surface': SurfaceFilter, 'universe': UniverseFilter, + 'zernike': ZernikeFilter } diff --git a/openmc/checkvalue.py b/openmc/checkvalue.py index 3b334319f6..cbeac9c373 100644 --- a/openmc/checkvalue.py +++ b/openmc/checkvalue.py @@ -265,8 +265,9 @@ def check_filetype_version(obj, expected_type, expected_version): if this_version[0] != expected_version: raise IOError('{} file has a version of {} which is not ' 'consistent with the version expected by OpenMC, {}' - .format(this_filetype, '.'.join(this_version), - expected_version)) + .format(this_filetype, + '.'.join(str(v) for v in this_version), + expected_version)) except AttributeError: raise IOError('Could not read {} file. This most likely means the {} ' 'file was produced by a different version of OpenMC than ' diff --git a/openmc/deplete/integrator/cecm.py b/openmc/deplete/integrator/cecm.py index 61b58d0b95..0b17b49b57 100644 --- a/openmc/deplete/integrator/cecm.py +++ b/openmc/deplete/integrator/cecm.py @@ -47,11 +47,42 @@ def cecm(operator, timesteps, power, print_out=True): # Generate initial conditions with operator as vec: chain = operator.chain - t = 0.0 + + # Initialize time + if operator.prev_res is None: + t = 0.0 + else: + t = operator.prev_res[-1].time[-1] + + # Initialize starting index for saving results + if operator.prev_res is None: + i_res = 0 + else: + i_res = len(operator.prev_res) + for i, (dt, p) in enumerate(zip(timesteps, power)): - # Get beginning-of-timestep reaction rates - x = [copy.deepcopy(vec)] - op_results = [operator(x[0], p)] + # Get beginning-of-timestep concentrations and reaction rates + # Avoid doing first transport run if already done in previous + # calculation + if i > 0 or operator.prev_res is None: + x = [copy.deepcopy(vec)] + op_results = [operator(x[0], p)] + + else: + # Get initial concentration + x = [operator.prev_res[-1].data[0]] + + # Get rates + op_results = [operator.prev_res[-1]] + op_results[0].rates = op_results[0].rates[0] + + # Set first stage value of keff + op_results[0].k = op_results[0].k[0] + + # Scale reaction rates by ratio of powers + power_res = operator.prev_res[-1].power + ratio_power = p / power_res + op_results[0].rates[0] *= ratio_power[0] # Deplete for first half of timestep x_middle = deplete(chain, x[0], op_results[0], dt/2, print_out) @@ -61,10 +92,11 @@ def cecm(operator, timesteps, power, print_out=True): op_results.append(operator(x_middle, p)) # Deplete for full timestep using beginning-of-step materials + # and middle-of-timestep reaction rates x_end = deplete(chain, x[0], op_results[1], dt, print_out) # Create results, write to disk - Results.save(operator, x, op_results, [t, t + dt], i) + Results.save(operator, x, op_results, [t, t + dt], p, i_res + i) # Advance time, update vector t += dt @@ -75,4 +107,4 @@ def cecm(operator, timesteps, power, print_out=True): op_results = [operator(x[0], power[-1])] # Create results, write to disk - Results.save(operator, x, op_results, [t, t], len(timesteps)) + Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps)) diff --git a/openmc/deplete/integrator/predictor.py b/openmc/deplete/integrator/predictor.py index 9be992c16a..a450178134 100644 --- a/openmc/deplete/integrator/predictor.py +++ b/openmc/deplete/integrator/predictor.py @@ -42,14 +42,41 @@ def predictor(operator, timesteps, power, print_out=True): # Generate initial conditions with operator as vec: chain = operator.chain - t = 0.0 - for i, (dt, p) in enumerate(zip(timesteps, power)): - # Get beginning-of-timestep reaction rates - x = [copy.deepcopy(vec)] - op_results = [operator(x[0], p)] - # Create results, write to disk - Results.save(operator, x, op_results, [t, t + dt], i) + # Initialize time + if operator.prev_res is None: + t = 0.0 + else: + t = operator.prev_res[-1].time[-1] + + # Initialize starting index for saving results + if operator.prev_res is None: + i_res = 0 + else: + i_res = len(operator.prev_res) - 1 + + for i, (dt, p) in enumerate(zip(timesteps, power)): + # Get beginning-of-timestep concentrations and reaction rates + # Avoid doing first transport run if already done in previous + # calculation + if i > 0 or operator.prev_res is None: + x = [copy.deepcopy(vec)] + op_results = [operator(x[0], p)] + + # Create results, write to disk + Results.save(operator, x, op_results, [t, t + dt], p, i_res + i) + else: + # Get initial concentration + x = [operator.prev_res[-1].data[0]] + + # Get rates + op_results = [operator.prev_res[-1]] + op_results[0].rates = op_results[0].rates[0] + + # Scale reaction rates by ratio of powers + power_res = operator.prev_res[-1].power + ratio_power = p / power_res + op_results[0].rates[0] *= ratio_power[0] # Deplete for full timestep x_end = deplete(chain, x[0], op_results[0], dt, print_out) @@ -63,4 +90,4 @@ def predictor(operator, timesteps, power, print_out=True): op_results = [operator(x[0], power[-1])] # Create results, write to disk - Results.save(operator, x, op_results, [t, t], len(timesteps)) + Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps)) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index b40ff63447..a0b2417c8e 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -66,6 +66,10 @@ class Operator(TransportOperator): chain_file : str, optional Path to the depletion chain XML file. Defaults to the :envvar:`OPENMC_DEPLETE_CHAIN` environment variable if it exists. + prev_results : ResultsList, optional + Results from a previous depletion calculation. If this argument is + specified, the depletion calculation will start from the latest state + in the previous results. Attributes ---------- @@ -75,7 +79,7 @@ class Operator(TransportOperator): OpenMC settings object dilute_initial : float Initial atom density to add for nuclides that are zero in initial - condition to ensure they exist in the decay chain. Only done for + condition to ensure they exist in the decay chain. Only done for nuclides with reaction rates. Defaults to 1.0e3. output_dir : pathlib.Path Path to output directory to save results. @@ -94,14 +98,25 @@ class Operator(TransportOperator): All burnable material IDs local_mats : list of str All burnable material IDs being managed by a single process + prev_res : ResultsList + Results from a previous depletion calculation """ - def __init__(self, geometry, settings, chain_file=None): + def __init__(self, geometry, settings, chain_file=None, prev_results=None): super().__init__(chain_file) self.round_number = False self.settings = settings self.geometry = geometry + if prev_results != None: + # Reload volumes into geometry + prev_results[-1].transfer_volumes(geometry) + + # Store previous results in operator + self.prev_res = prev_results + else: + self.prev_res = None + # Clear out OpenMC, create task lists, distribute openmc.reset_auto_ids() self.burnable_mats, volume, nuclides = self._get_burnable_mats() @@ -109,16 +124,19 @@ class Operator(TransportOperator): # Determine which nuclides have incident neutron data self.nuclides_with_data = self._get_nuclides_with_data() - self._burnable_nucs = [nuc for nuc in self.nuclides_with_data - if nuc in self.chain] - # Extract number densities from the geometry - self._extract_number(self.local_mats, volume, nuclides) + # Select nuclides with data that are also in the chain + self._burnable_nucs = [nuc.name for nuc in self.chain.nuclides + if nuc.name in self.nuclides_with_data] + + # Extract number densities from the geometry / previous depletion run + self._extract_number(self.local_mats, volume, nuclides, self.prev_res) # Create reaction rates array self.reaction_rates = ReactionRates( self.local_mats, self._burnable_nucs, self.chain.reactions) + def __call__(self, vec, power, print_out=True): """Runs a simulation. @@ -212,7 +230,7 @@ class Operator(TransportOperator): return burnable_mats, volume, nuclides - def _extract_number(self, local_mats, volume, nuclides): + def _extract_number(self, local_mats, volume, nuclides, prev_res=None): """Construct AtomNumber using geometry Parameters @@ -223,6 +241,8 @@ class Operator(TransportOperator): Volumes for the above materials in [cm^3] nuclides : list of str Nuclides to be used in the simulation. + prev_res : ResultsList, optional + Results from a previous depletion calculation """ self.number = AtomNumber(local_mats, nuclides, volume, len(self.chain)) @@ -231,10 +251,18 @@ class Operator(TransportOperator): for nuc in self._burnable_nucs: self.number.set_atom_density(np.s_[:], nuc, self.dilute_initial) - # Now extract the number densities and store - for mat in self.geometry.get_all_materials().values(): - if str(mat.id) in local_mats: - self._set_number_from_mat(mat) + # Now extract and store the number densities + # From the geometry if no previous depletion results + if prev_res is None: + for mat in self.geometry.get_all_materials().values(): + if str(mat.id) in local_mats: + self._set_number_from_mat(mat) + + # Else from previous depletion results + else: + for mat in self.geometry.get_all_materials().values(): + if str(mat.id) in local_mats: + self._set_number_from_results(mat, prev_res) def _set_number_from_mat(self, mat): """Extracts material and number densities from openmc.Material @@ -251,6 +279,41 @@ class Operator(TransportOperator): number = density * 1.0e24 self.number.set_atom_density(mat_id, nuclide, number) + def _set_number_from_results(self, mat, prev_res): + """Extracts material nuclides and number densities. + + If the nuclide concentration's evolution is tracked, the densities come + from depletion results. Else, densities are extracted from the geometry + in the summary. + + Parameters + ---------- + mat : openmc.Material + The material to read from + prev_res : ResultsList + Results from a previous depletion calculation + + """ + mat_id = str(mat.id) + + # Get nuclide lists from geometry and depletion results + depl_nuc = prev_res[-1].nuc_to_ind + geom_nuc_densities = mat.get_nuclide_atom_densities() + + # Merge lists of nuclides, with the same order for every calculation + geom_nuc_densities.update(depl_nuc) + + for nuclide in geom_nuc_densities.keys(): + if nuclide in depl_nuc: + concentration = prev_res.get_atoms(mat_id, nuclide)[1][-1] + volume = prev_res[-1].volume[mat_id] + number = concentration / volume + else: + density = geom_nuc_densities[nuclide][1] + number = density * 1.0e24 + + self.number.set_atom_density(mat_id, nuclide, number) + def initial_condition(self): """Performs final setup and returns initial condition. @@ -306,15 +369,19 @@ class Operator(TransportOperator): densities.append(val) else: # Only output warnings if values are significantly - # negative. CRAM does not guarantee positive values. + # negative. CRAM does not guarantee positive values. if val < -1.0e-21: print("WARNING: nuclide ", nuc, " in material ", mat, " is negative (density = ", val, " at/barn-cm)") number_i[mat, nuc] = 0.0 + # Update densities on C API side mat_internal = openmc.capi.materials[int(mat)] mat_internal.set_densities(nuclides, densities) + #TODO Update densities on the Python side, otherwise the + # summary.h5 file contains densities at the first time step + def _generate_materials_xml(self): """Creates materials.xml from self.number. diff --git a/openmc/deplete/reaction_rates.py b/openmc/deplete/reaction_rates.py index fddb88b19d..cea2f19976 100644 --- a/openmc/deplete/reaction_rates.py +++ b/openmc/deplete/reaction_rates.py @@ -22,6 +22,9 @@ class ReactionRates(np.ndarray): Depletable nuclides reactions : list of str Transmutation reactions being tracked + from_results : boolean + If the reaction rates are loaded from results, indexing dictionnaries + need to be kept the same. Attributes ---------- @@ -47,16 +50,24 @@ class ReactionRates(np.ndarray): # the __array_finalize__ method (discussed here: # https://docs.scipy.org/doc/numpy/user/basics.subclassing.html) - def __new__(cls, local_mats, nuclides, reactions): + def __new__(cls, local_mats, nuclides, reactions, from_results=False): # Create appropriately-sized zeroed-out ndarray shape = (len(local_mats), len(nuclides), len(reactions)) obj = super().__new__(cls, shape) obj[:] = 0.0 - # Add mapping attributes - obj.index_mat = {mat: i for i, mat in enumerate(local_mats)} - obj.index_nuc = {nuc: i for i, nuc in enumerate(nuclides)} - obj.index_rx = {rx: i for i, rx in enumerate(reactions)} + # Add mapping attributes, keep same indexing if from depletion_results + if from_results: + obj.index_mat = local_mats + obj.index_nuc = nuclides + obj.index_rx = reactions + # Else, assumes that reaction rates are ordered the same way as + # the lists of local_mats, nuclides and reactions (or keys if these + # are dictionnaries) + else: + obj.index_mat = {mat: i for i, mat in enumerate(local_mats)} + obj.index_nuc = {nuc: i for i, nuc in enumerate(nuclides)} + obj.index_rx = {rx: i for i, rx in enumerate(reactions)} return obj diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py index ab740e61ec..6170a44643 100644 --- a/openmc/deplete/results.py +++ b/openmc/deplete/results.py @@ -5,6 +5,7 @@ Contains results generation and saving capabilities. from collections import OrderedDict import copy +from warnings import warn import numpy as np import h5py @@ -24,6 +25,8 @@ class Results(object): Eigenvalue for each substep. time : list of float Time at beginning, end of step, in seconds. + power : float + Power during time step, in Watts n_mat : int Number of mats. n_nuc : int @@ -49,6 +52,7 @@ class Results(object): def __init__(self): self.k = None self.time = None + self.power = None self.rates = None self.volume = None @@ -161,6 +165,7 @@ class Results(object): else: kwargs = {} + # Write new file if first time step, else add to existing file kwargs['mode'] = "w" if step == 0 else "a" with h5py.File(filename, **kwargs) as handle: @@ -237,6 +242,9 @@ class Results(object): handle.create_dataset("time", (1, 2), maxshape=(None, 2), dtype='float64') + handle.create_dataset("power", (1, n_stages), maxshape=(None, n_stages), + dtype='float64') + def _to_hdf5(self, handle, index): """Converts results object into an hdf5 object. @@ -259,6 +267,7 @@ class Results(object): rxn_dset = handle["/reaction rates"] eigenvalues_dset = handle["/eigenvalues"] time_dset = handle["/time"] + power_dset = handle["/power"] # Get number of results stored number_shape = list(number_dset.shape) @@ -283,6 +292,10 @@ class Results(object): time_shape[0] = new_shape time_dset.resize(time_shape) + power_shape = list(power_dset.shape) + power_shape[0] = new_shape + power_dset.resize(power_shape) + # If nothing to write, just return if len(self.mat_to_ind) == 0: return @@ -300,6 +313,7 @@ class Results(object): eigenvalues_dset[index, i] = self.k[i] if comm.rank == 0: time_dset[index, :] = self.time + power_dset[index, :] = self.power @classmethod def from_hdf5(cls, handle, step): @@ -319,10 +333,12 @@ class Results(object): number_dset = handle["/number"] eigenvalues_dset = handle["/eigenvalues"] time_dset = handle["/time"] + power_dset = handle["/power"] results.data = number_dset[step, :, :, :] results.k = eigenvalues_dset[step, :] results.time = time_dset[step, :] + results.power = power_dset[step, :] # Reconstruct dictionaries results.volume = OrderedDict() @@ -351,7 +367,7 @@ class Results(object): results.rates = [] # Reconstruct reactions for i in range(results.n_stages): - rate = ReactionRates(results.mat_to_ind, rxn_nuc_to_ind, rxn_to_ind) + rate = ReactionRates(results.mat_to_ind, rxn_nuc_to_ind, rxn_to_ind, True) rate[:] = handle["/reaction rates"][step, i, :, :, :] results.rates.append(rate) @@ -359,7 +375,7 @@ class Results(object): return results @staticmethod - def save(op, x, op_results, t, step_ind): + def save(op, x, op_results, t, power, step_ind): """Creates and writes depletion results to disk Parameters @@ -372,6 +388,8 @@ class Results(object): Results of applying transport operator t : list of float Time indices. + power : float + Power during time step step_ind : int Step index. @@ -379,8 +397,18 @@ class Results(object): # Get indexing terms vol_dict, nuc_list, burn_list, full_burn_list = op.get_results_info() - # Create results + # For a restart calculation, limit number of stages saved to meet the + # format of the hdf5 file stages = len(x) + offset = 0 + if op.prev_res is not None and op.prev_res[0].n_stages < stages: + offset = stages - op.prev_res[0].n_stages + stages = min(stages, op.prev_res[0].n_stages) + warn("Number of restart integrator stages saved limited by initial" + " depletion integrator choice to {}" + .format(op.prev_res[0].n_stages)) + + # Create results results = Results() results.allocate(vol_dict, nuc_list, burn_list, full_burn_list, stages) @@ -388,10 +416,25 @@ class Results(object): for i in range(stages): for mat_i in range(n_mat): - results[i, mat_i, :] = x[i][mat_i][:] + results[i, mat_i, :] = x[offset + i][mat_i][:] results.k = [r.k for r in op_results] results.rates = [r.rates for r in op_results] results.time = t + results.power = power results.export_to_hdf5("depletion_results.h5", step_ind) + + def transfer_volumes(self, geometry): + """Transfers volumes from depletion results to geometry + + Parameters + ---------- + geometry : OpenMC geometry to be used in a depletion restart + calculation + + """ + for cell in geometry.get_all_material_cells().values(): + for material in cell.get_all_materials().values(): + if material.depletable: + material.volume = self.volume[str(material.id)] diff --git a/openmc/deplete/results_list.py b/openmc/deplete/results_list.py index d3d45e955b..58e4b66dc1 100644 --- a/openmc/deplete/results_list.py +++ b/openmc/deplete/results_list.py @@ -43,8 +43,8 @@ class ResultsList(list): Total number of atoms for specified nuclide """ - time = np.empty_like(self) - concentration = np.empty_like(self) + time = np.empty_like(self, dtype=float) + concentration = np.empty_like(self, dtype=float) # Evaluate value in each region for i, result in enumerate(self): @@ -73,8 +73,8 @@ class ResultsList(list): Array of reaction rates """ - time = np.empty_like(self) - rate = np.empty_like(self) + time = np.empty_like(self, dtype=float) + rate = np.empty_like(self, dtype=float) # Evaluate value in each region for i, result in enumerate(self): @@ -94,8 +94,8 @@ class ResultsList(list): k-eigenvalue at each time """ - time = np.empty_like(self) - eigenvalue = np.empty_like(self) + time = np.empty_like(self, dtype=float) + eigenvalue = np.empty_like(self, dtype=float) # Get time/eigenvalue at each point for i, result in enumerate(self): diff --git a/openmc/lattice.py b/openmc/lattice.py index 86cfd5c41b..ef33247fca 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -546,10 +546,15 @@ class RectLattice(Lattice): """ if self.ndim == 2: nx, ny = self.shape - return np.broadcast(*np.ogrid[:nx, :ny]) + for iy in range(ny): + for ix in range(nx): + yield (ix, iy) else: nx, ny, nz = self.shape - return np.broadcast(*np.ogrid[:nx, :ny, :nz]) + for iz in range(nz): + for iy in range(ny): + for ix in range(nx): + yield (ix, iy, iz) @property def lower_left(self): diff --git a/openmc/material.py b/openmc/material.py index cfa3eba3c3..dac2fe14af 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -467,7 +467,7 @@ class Material(IDManagerMixin): raise ValueError(msg) # Generally speaking, the density for a macroscopic object will - # be 1.0. Therefore, lets set density to 1.0 so that the user + # be 1.0. Therefore, lets set density to 1.0 so that the user # doesnt need to set it unless its needed. # Of course, if the user has already set a value of density, # then we will not override it. diff --git a/openmc/mesh.py b/openmc/mesh.py index c9c76552bc..3f53582d5e 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -182,6 +182,40 @@ class Mesh(IDManagerMixin): return mesh + @classmethod + def from_rect_lattice(cls, lattice, division=1, mesh_id=None, name=''): + """Create mesh from an existing rectangular lattice + + Parameters + ---------- + lattice : openmc.RectLattice + Rectangular lattice used as a template for this mesh + division : int + Number of mesh cells per lattice cell. + If not specified, there will be 1 mesh cell per lattice cell. + mesh_id : int + Unique identifier for the mesh + name : str + Name of the mesh + + Returns + ------- + openmc.Mesh + Mesh instance + + """ + cv.check_type('rectangular lattice', lattice, openmc.RectLattice) + + shape = np.array(lattice.shape) + width = lattice.pitch*shape + + mesh = cls(mesh_id, name) + mesh.lower_left = lattice.lower_left + mesh.upper_right = lattice.lower_left + width + mesh.dimension = shape*division + + return mesh + def to_xml_element(self): """Return XML representation of the mesh diff --git a/openmc/tallies.py b/openmc/tallies.py index a5341cbedd..74f3429034 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -248,7 +248,7 @@ class Tally(IDManagerMixin): @property def sum_sq(self): - if not self._sp_filename: + if not self._sp_filename or self.derived: return None if not self._results_read: diff --git a/src/api.F90 b/src/api.F90 index f32de6bebb..ec17472ea2 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -35,6 +35,7 @@ module openmc_api private public :: openmc_calculate_volumes + public :: openmc_cell_filter_get_bins public :: openmc_cell_get_id public :: openmc_cell_get_fill public :: openmc_cell_set_fill @@ -206,7 +207,7 @@ contains if (found) then if (rtype == 1) then - id = cells(p % coord(p % n_coord) % cell) % id + id = cells(p % coord(p % n_coord) % cell) % id() elseif (rtype == 2) then if (p % material == MATERIAL_VOID) then id = 0 diff --git a/src/cell.cpp b/src/cell.cpp new file mode 100644 index 0000000000..e4388e7012 --- /dev/null +++ b/src/cell.cpp @@ -0,0 +1,523 @@ +#include "cell.h" + +#include +#include +#include +#include + +#include "constants.h" +#include "error.h" +#include "hdf5_interface.h" +#include "lattice.h" +#include "surface.h" +#include "xml_interface.h" + + +namespace openmc { + +//============================================================================== +// Constants +//============================================================================== + +constexpr int32_t OP_LEFT_PAREN {std::numeric_limits::max()}; +constexpr int32_t OP_RIGHT_PAREN {std::numeric_limits::max() - 1}; +constexpr int32_t OP_COMPLEMENT {std::numeric_limits::max() - 2}; +constexpr int32_t OP_INTERSECTION {std::numeric_limits::max() - 3}; +constexpr int32_t OP_UNION {std::numeric_limits::max() - 4}; + +extern "C" double FP_PRECISION; + +//============================================================================== +// Global variables +//============================================================================== + +int32_t n_cells {0}; + +std::vector global_cells; +std::unordered_map cell_map; + +std::vector global_universes; +std::unordered_map universe_map; + +//============================================================================== +//! Convert region specification string to integer tokens. +//! +//! The characters (, ), |, and ~ count as separate tokens since they represent +//! operators. +//============================================================================== + +std::vector +tokenize(const std::string region_spec) { + // Check for an empty region_spec first. + std::vector tokens; + if (region_spec.empty()) { + return tokens; + } + + // Parse all halfspaces and operators except for intersection (whitespace). + for (int i = 0; i < region_spec.size(); ) { + if (region_spec[i] == '(') { + tokens.push_back(OP_LEFT_PAREN); + i++; + + } else if (region_spec[i] == ')') { + tokens.push_back(OP_RIGHT_PAREN); + i++; + + } else if (region_spec[i] == '|') { + tokens.push_back(OP_UNION); + i++; + + } else if (region_spec[i] == '~') { + tokens.push_back(OP_COMPLEMENT); + i++; + + } else if (region_spec[i] == '-' || region_spec[i] == '+' + || std::isdigit(region_spec[i])) { + // This is the start of a halfspace specification. Iterate j until we + // find the end, then push-back everything between i and j. + int j = i + 1; + while (j < region_spec.size() && std::isdigit(region_spec[j])) {j++;} + tokens.push_back(std::stoi(region_spec.substr(i, j-i))); + i = j; + + } else if (std::isspace(region_spec[i])) { + i++; + + } else { + std::stringstream err_msg; + err_msg << "Region specification contains invalid character, \"" + << region_spec[i] << "\""; + fatal_error(err_msg); + } + } + + // Add in intersection operators where a missing operator is needed. + int i = 0; + while (i < tokens.size()-1) { + bool left_compat {(tokens[i] < OP_UNION) || (tokens[i] == OP_RIGHT_PAREN)}; + bool right_compat {(tokens[i+1] < OP_UNION) + || (tokens[i+1] == OP_LEFT_PAREN) + || (tokens[i+1] == OP_COMPLEMENT)}; + if (left_compat && right_compat) { + tokens.insert(tokens.begin()+i+1, OP_INTERSECTION); + } + i++; + } + + return tokens; +} + +//============================================================================== +//! Convert infix region specification to Reverse Polish Notation (RPN) +//! +//! This function uses the shunting-yard algorithm. +//============================================================================== + +std::vector +generate_rpn(int32_t cell_id, std::vector infix) +{ + std::vector rpn; + std::vector stack; + + for (int32_t token : infix) { + if (token < OP_UNION) { + // If token is not an operator, add it to output + rpn.push_back(token); + + } else if (token < OP_RIGHT_PAREN) { + // Regular operators union, intersection, complement + while (stack.size() > 0) { + int32_t op = stack.back(); + + if (op < OP_RIGHT_PAREN && + ((token == OP_COMPLEMENT && token < op) || + (token != OP_COMPLEMENT && token <= op))) { + // While there is an operator, op, on top of the stack, if the token + // is left-associative and its precedence is less than or equal to + // that of op or if the token is right-associative and its precedence + // is less than that of op, move op to the output queue and push the + // token on to the stack. Note that only complement is + // right-associative. + rpn.push_back(op); + stack.pop_back(); + } else { + break; + } + } + + stack.push_back(token); + + } else if (token == OP_LEFT_PAREN) { + // If the token is a left parenthesis, push it onto the stack + stack.push_back(token); + + } else { + // If the token is a right parenthesis, move operators from the stack to + // the output queue until reaching the left parenthesis. + for (auto it = stack.rbegin(); *it != OP_LEFT_PAREN; it++) { + // If we run out of operators without finding a left parenthesis, it + // means there are mismatched parentheses. + if (it == stack.rend()) { + std::stringstream err_msg; + err_msg << "Mismatched parentheses in region specification for cell " + << cell_id; + fatal_error(err_msg); + } + + rpn.push_back(stack.back()); + stack.pop_back(); + } + + // Pop the left parenthesis. + stack.pop_back(); + } + } + + while (stack.size() > 0) { + int32_t op = stack.back(); + + // If the operator is a parenthesis it is mismatched. + if (op >= OP_RIGHT_PAREN) { + std::stringstream err_msg; + err_msg << "Mismatched parentheses in region specification for cell " + << cell_id; + fatal_error(err_msg); + } + + rpn.push_back(stack.back()); + stack.pop_back(); + } + + return rpn; +} + +//============================================================================== +// Cell implementation +//============================================================================== + +Cell::Cell(pugi::xml_node cell_node) +{ + if (check_for_node(cell_node, "id")) { + id = stoi(get_node_value(cell_node, "id")); + } else { + fatal_error("Must specify id of cell in geometry XML file."); + } + + //TODO: don't automatically lowercase cell and surface names + if (check_for_node(cell_node, "name")) { + name = get_node_value(cell_node, "name"); + } + + if (check_for_node(cell_node, "universe")) { + universe = stoi(get_node_value(cell_node, "universe")); + } else { + universe = 0; + } + + if (check_for_node(cell_node, "fill")) { + fill = stoi(get_node_value(cell_node, "fill")); + } else { + fill = C_NONE; + } + + if (check_for_node(cell_node, "material")) { + //TODO: read material ids. + material.push_back(C_NONE+1); + material.shrink_to_fit(); + } else { + material.push_back(C_NONE); + material.shrink_to_fit(); + } + + // Make sure that either material or fill was specified. + if ((material[0] == C_NONE) && (fill == C_NONE)) { + std::stringstream err_msg; + err_msg << "Neither material nor fill was specified for cell " << id; + fatal_error(err_msg); + } + + // Make sure that material and fill haven't been specified simultaneously. + if ((material[0] != C_NONE) && (fill != C_NONE)) { + std::stringstream err_msg; + err_msg << "Cell " << id << " has both a material and a fill specified; " + << "only one can be specified per cell"; + fatal_error(err_msg); + } + + // Read the region specification. + std::string region_spec; + if (check_for_node(cell_node, "region")) { + region_spec = get_node_value(cell_node, "region"); + } + + // Get a tokenized representation of the region specification. + region = tokenize(region_spec); + region.shrink_to_fit(); + + // Convert user IDs to surface indices. + for (auto &r : region) { + if (r < OP_UNION) { + r = copysign(surface_map[abs(r)] + 1, r); + } + } + + // Convert the infix region spec to RPN. + rpn = generate_rpn(id, region); + rpn.shrink_to_fit(); + + // Check if this is a simple cell. + simple = true; + for (int32_t token : rpn) { + if ((token == OP_COMPLEMENT) || (token == OP_UNION)) { + simple = false; + break; + } + } +} + +//============================================================================== + +bool +Cell::contains(const double xyz[3], const double uvw[3], + int32_t on_surface) const +{ + if (simple) { + return contains_simple(xyz, uvw, on_surface); + } else { + return contains_complex(xyz, uvw, on_surface); + } +} + +//============================================================================== + +std::pair +Cell::distance(const double xyz[3], const double uvw[3], + int32_t on_surface) const +{ + double min_dist {INFTY}; + int32_t i_surf {std::numeric_limits::max()}; + + for (int32_t token : rpn) { + // Ignore this token if it corresponds to an operator rather than a region. + if (token >= OP_UNION) {continue;} + + // Calculate the distance to this surface. + // Note the off-by-one indexing + bool coincident {token == on_surface}; + double d {surfaces_c[abs(token)-1]->distance(xyz, uvw, coincident)}; + + // Check if this distance is the new minimum. + if (d < min_dist) { + if (std::abs(d - min_dist) / min_dist >= FP_PRECISION) { + min_dist = d; + i_surf = -token; + } + } + } + + return {min_dist, i_surf}; +} + +//============================================================================== + +void +Cell::to_hdf5(hid_t cell_group) const +{ + if (!name.empty()) { + write_string(cell_group, "name", name, false); + } + + //TODO: Fix the off-by-one indexing. + write_int(cell_group, 0, nullptr, "universe", + &global_universes[universe-1]->id, false); + + // Write the region specification. + if (!region.empty()) { + std::stringstream region_spec {}; + for (int32_t token : region) { + if (token == OP_LEFT_PAREN) { + region_spec << " ("; + } else if (token == OP_RIGHT_PAREN) { + region_spec << " )"; + } else if (token == OP_COMPLEMENT) { + region_spec << " ~"; + } else if (token == OP_INTERSECTION) { + } else if (token == OP_UNION) { + region_spec << " |"; + } else { + // Note the off-by-one indexing + region_spec << " " << copysign(surfaces_c[abs(token)-1]->id, token); + } + } + write_string(cell_group, "region", region_spec.str(), false); + } +} + +//============================================================================== + +bool +Cell::contains_simple(const double xyz[3], const double uvw[3], + int32_t on_surface) const +{ + for (int32_t token : rpn) { + if (token < OP_UNION) { + // If the token is not an operator, evaluate the sense of particle with + // respect to the surface and see if the token matches the sense. If the + // particle's surface attribute is set and matches the token, that + // overrides the determination based on sense(). + if (token == on_surface) { + } else if (-token == on_surface) { + return false; + } else { + // Note the off-by-one indexing + bool sense = surfaces_c[abs(token)-1]->sense(xyz, uvw); + if (sense != (token > 0)) {return false;} + } + } + } + return true; +} + +//============================================================================== + +bool +Cell::contains_complex(const double xyz[3], const double uvw[3], + int32_t on_surface) const +{ + // Make a stack of booleans. We don't know how big it needs to be, but we do + // know that rpn.size() is an upper-bound. + bool stack[rpn.size()]; + int i_stack = -1; + + for (int32_t token : rpn) { + // If the token is a binary operator (intersection/union), apply it to + // the last two items on the stack. If the token is a unary operator + // (complement), apply it to the last item on the stack. + if (token == OP_UNION) { + stack[i_stack-1] = stack[i_stack-1] || stack[i_stack]; + i_stack --; + } else if (token == OP_INTERSECTION) { + stack[i_stack-1] = stack[i_stack-1] && stack[i_stack]; + i_stack --; + } else if (token == OP_COMPLEMENT) { + stack[i_stack] = !stack[i_stack]; + } else { + // If the token is not an operator, evaluate the sense of particle with + // respect to the surface and see if the token matches the sense. If the + // particle's surface attribute is set and matches the token, that + // overrides the determination based on sense(). + i_stack ++; + if (token == on_surface) { + stack[i_stack] = true; + } else if (-token == on_surface) { + stack[i_stack] = false; + } else { + // Note the off-by-one indexing + bool sense = surfaces_c[abs(token)-1]->sense(xyz, uvw);; + stack[i_stack] = (sense == (token > 0)); + } + } + } + + if (i_stack == 0) { + // The one remaining bool on the stack indicates whether the particle is + // in the cell. + return stack[i_stack]; + } else { + // This case occurs if there is no region specification since i_stack will + // still be -1. + return true; + } +} + +//============================================================================== +// Non-method functions +//============================================================================== + +extern "C" void +read_cells(pugi::xml_node *node) +{ + // Count the number of cells. + for (pugi::xml_node cell_node: node->children("cell")) {n_cells++;} + if (n_cells == 0) { + fatal_error("No cells found in geometry.xml!"); + } + + // Allocate the vector of Cells. + global_cells.reserve(n_cells); + + // Loop over XML cell elements and populate the array. + for (pugi::xml_node cell_node: node->children("cell")) { + global_cells.push_back(new Cell(cell_node)); + } + + // Populate the Universe vector and map. + for (int i = 0; i < global_cells.size(); i++) { + int32_t uid = global_cells[i]->universe; + auto it = universe_map.find(uid); + if (it == universe_map.end()) { + global_universes.push_back(new Universe()); + global_universes.back()->id = uid; + global_universes.back()->cells.push_back(i); + universe_map[uid] = global_universes.size() - 1; + } else { + global_universes[it->second]->cells.push_back(i); + } + } +} + +//============================================================================== +// Fortran compatibility functions +//============================================================================== + +extern "C" { + Cell* cell_pointer(int32_t cell_ind) {return global_cells[cell_ind];} + + int32_t cell_id(Cell *c) {return c->id;} + + void cell_set_id(Cell *c, int32_t id) {c->id = id;} + + int cell_type(Cell *c) {return c->type;} + + void cell_set_type(Cell *c, int type) {c->type = type;} + + int32_t cell_universe(Cell *c) {return c->universe;} + + void cell_set_universe(Cell *c, int32_t universe) {c->universe = universe;} + + int32_t cell_fill(Cell *c) {return c->fill;} + + int32_t* cell_fill_ptr(Cell *c) {return &c->fill;} + + int32_t cell_n_instances(Cell *c) {return c->n_instances;} + + bool cell_simple(Cell *c) {return c->simple;} + + bool cell_contains(Cell *c, double xyz[3], double uvw[3], int32_t on_surface) + {return c->contains(xyz, uvw, on_surface);} + + void cell_distance(Cell *c, double xyz[3], double uvw[3], int32_t on_surface, + double *min_dist, int32_t *i_surf) + { + std::pair out = c->distance(xyz, uvw, on_surface); + *min_dist = out.first; + *i_surf = out.second; + } + + int32_t cell_offset(Cell *c, int map) {return c->offset[map];} + + void cell_to_hdf5(Cell *c, hid_t group) {c->to_hdf5(group);} + + void extend_cells_c(int32_t n) + { + global_cells.reserve(global_cells.size() + n); + for (int32_t i = 0; i < n; i++) { + global_cells.push_back(new Cell()); + } + n_cells = global_cells.size(); + } +} + + +} // namespace openmc diff --git a/src/cell.h b/src/cell.h new file mode 100644 index 0000000000..a2f72c8337 --- /dev/null +++ b/src/cell.h @@ -0,0 +1,118 @@ +#ifndef CELL_H +#define CELL_H + +#include +#include +#include +#include + +#include "hdf5.h" +#include "pugixml/pugixml.hpp" + + +namespace openmc { + +//============================================================================== +// Constants +//============================================================================== + +extern "C" int FILL_MATERIAL; +extern "C" int FILL_UNIVERSE; +extern "C" int FILL_LATTICE; + +//============================================================================== +// Global variables +//============================================================================== + +extern "C" int32_t n_cells; + +class Cell; +extern std::vector global_cells; +extern std::unordered_map cell_map; + +class Universe; +extern std::vector global_universes; +extern std::unordered_map universe_map; + +//============================================================================== +//! A geometry primitive that fills all space and contains cells. +//============================================================================== + +class Universe +{ +public: + int32_t id; //!< Unique ID + std::vector cells; //!< Cells within this universe + //double x0, y0, z0; //!< Translation coordinates. +}; + +//============================================================================== +//! A geometry primitive that links surfaces, universes, and materials +//============================================================================== + +class Cell +{ +public: + int32_t id; //!< Unique ID + std::string name; //!< User-defined name + int type; //!< Material, universe, or lattice + int32_t universe; //!< Universe # this cell is in + int32_t fill; //!< Universe # filling this cell + int32_t n_instances{0}; //!< Number of instances of this cell + + //! \brief Material(s) within this cell. + //! + //! May be multiple materials for distribcell. C_NONE signifies a universe. + std::vector material; + + //! Definition of spatial region as Boolean expression of half-spaces + std::vector region; + //! Reverse Polish notation for region expression + std::vector rpn; + bool simple; //!< Does the region contain only intersections? + + std::vector offset; //!< Distribcell offset table + + Cell() {}; + + explicit Cell(pugi::xml_node cell_node); + + //! \brief Determine if a cell contains the particle at a given location. + //! + //! The bounds of the cell are detemined by a logical expression involving + //! surface half-spaces. At initialization, the expression was converted + //! to RPN notation. + //! + //! The function is split into two cases, one for simple cells (those + //! involving only the intersection of half-spaces) and one for complex cells. + //! Simple cells can be evaluated with short circuit evaluation, i.e., as soon + //! as we know that one half-space is not satisfied, we can exit. This + //! provides a performance benefit for the common case. In + //! contains_complex, we evaluate the RPN expression using a stack, similar to + //! how a RPN calculator would work. + //! @param xyz[3] The 3D Cartesian coordinate to check. + //! @param uvw[3] A direction used to "break ties" the coordinates are very + //! close to a surface. + //! @param on_surface The signed index of a surface that the coordinate is + //! known to be on. This index takes precedence over surface sense + //! calculations. + bool + contains(const double xyz[3], const double uvw[3], int32_t on_surface) const; + + //! Find the oncoming boundary of this cell. + std::pair + distance(const double xyz[3], const double uvw[3], int32_t on_surface) const; + + //! \brief Write cell information to an HDF5 group. + //! @param group_id An HDF5 group id. + void to_hdf5(hid_t group_id) const; + +protected: + bool contains_simple(const double xyz[3], const double uvw[3], + int32_t on_surface) const; + bool contains_complex(const double xyz[3], const double uvw[3], + int32_t on_surface) const; +}; + +} // namespace openmc +#endif // CELL_H diff --git a/src/constants.F90 b/src/constants.F90 index 430cc27151..035494b87b 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -113,6 +113,9 @@ module constants FILL_MATERIAL = 1, & ! Cell with a specified material FILL_UNIVERSE = 2, & ! Cell filled by a separate universe FILL_LATTICE = 3 ! Cell filled with a lattice + integer(C_INT), bind(C, name='FILL_MATERIAL') :: FILL_MATERIAL_C = FILL_MATERIAL + integer(C_INT), bind(C, name='FILL_UNIVERSE') :: FILL_UNIVERSE_C = FILL_UNIVERSE + integer(C_INT), bind(C, name='FILL_LATTICE') :: FILL_LATTICE_C = FILL_LATTICE ! Void material integer, parameter :: MATERIAL_VOID = -1 @@ -437,6 +440,7 @@ module constants ! indicates that an array index hasn't been set integer, parameter :: NONE = 0 + integer, parameter :: C_NONE = -1 ! Codes for read errors -- better hope these numbers are never used in an ! input file! diff --git a/src/constants.h b/src/constants.h index 3ddadb5a10..3b20e68f6a 100644 --- a/src/constants.h +++ b/src/constants.h @@ -6,6 +6,7 @@ #include #include +#include #include namespace openmc { @@ -80,7 +81,11 @@ constexpr int MG_GET_XS_NU_FISSION {12}; constexpr int MG_GET_XS_CHI_PROMPT {13}; constexpr int MG_GET_XS_CHI_DELAYED {14}; +extern "C" double FP_COINCIDENT; +extern "C" double FP_PRECISION; +constexpr double INFTY {std::numeric_limits::max()}; +constexpr int C_NONE {-1}; } // namespace openmc -#endif // CONSTANTS_H \ No newline at end of file +#endif // CONSTANTS_H diff --git a/src/geometry.F90 b/src/geometry.F90 index e747b619ac..457fdf3786 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -14,125 +14,36 @@ module geometry implicit none + interface + function cell_contains_c(cell_ptr, xyz, uvw, on_surface) & + bind(C, name="cell_contains") result(in_cell) + import C_PTR, C_DOUBLE, C_INT32_T, C_BOOL + type(C_PTR), intent(in), value :: cell_ptr + real(C_DOUBLE), intent(in) :: xyz(3) + real(C_DOUBLE), intent(in) :: uvw(3) + integer(C_INT32_T), intent(in), value :: on_surface + logical(C_BOOL) :: in_cell + end function cell_contains_c + + function count_universe_instances(search_univ, target_univ_id) bind(C) & + result(count) + import C_INT32_T, C_INT + integer(C_INT32_T), intent(in), value :: search_univ + integer(C_INT32_T), intent(in), value :: target_univ_id + integer(C_INT) :: count + end function + end interface + contains -!=============================================================================== -! CELL_CONTAINS determines if a cell contains the particle at a given -! location. The bounds of the cell are detemined by a logical expression -! involving surface half-spaces. At initialization, the expression was converted -! to RPN notation. -! -! The function is split into two cases, one for simple cells (those involving -! only the intersection of half-spaces) and one for complex cells. Simple cells -! can be evaluated with short circuit evaluation, i.e., as soon as we know that -! one half-space is not satisfied, we can exit. This provides a performance -! benefit for the common case. In complex_cell_contains, we evaluate the RPN -! expression using a stack, similar to how a RPN calculator would work. -!=============================================================================== - - pure function cell_contains(c, p) result(in_cell) + function cell_contains(c, p) result(in_cell) type(Cell), intent(in) :: c type(Particle), intent(in) :: p logical :: in_cell - - if (c%simple) then - in_cell = simple_cell_contains(c, p) - else - in_cell = complex_cell_contains(c, p) - end if + in_cell = cell_contains_c(c%ptr, p%coord(p%n_coord)%xyz, & + p%coord(p%n_coord)%uvw, p%surface) end function cell_contains - pure function simple_cell_contains(c, p) result(in_cell) - type(Cell), intent(in) :: c - type(Particle), intent(in) :: p - logical :: in_cell - - integer :: i - integer :: token - logical :: actual_sense ! sense of particle wrt surface - - in_cell = .true. - do i = 1, size(c%rpn) - token = c%rpn(i) - if (token < OP_UNION) then - ! If the token is not an operator, evaluate the sense of particle with - ! respect to the surface and see if the token matches the sense. If the - ! particle's surface attribute is set and matches the token, that - ! overrides the determination based on sense(). - if (token == p%surface) then - cycle - elseif (-token == p%surface) then - in_cell = .false. - exit - else - actual_sense = surfaces(abs(token)) % sense( & - p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw) - if (actual_sense .neqv. (token > 0)) then - in_cell = .false. - exit - end if - end if - end if - end do - end function simple_cell_contains - - pure function complex_cell_contains(c, p) result(in_cell) - type(Cell), intent(in) :: c - type(Particle), intent(in) :: p - logical :: in_cell - - integer :: i - integer :: token - integer :: i_stack - logical :: actual_sense ! sense of particle wrt surface - logical :: stack(size(c%rpn)) - - i_stack = 0 - do i = 1, size(c%rpn) - token = c%rpn(i) - - ! If the token is a binary operator (intersection/union), apply it to - ! the last two items on the stack. If the token is a unary operator - ! (complement), apply it to the last item on the stack. - select case (token) - case (OP_UNION) - stack(i_stack - 1) = stack(i_stack - 1) .or. stack(i_stack) - i_stack = i_stack - 1 - case (OP_INTERSECTION) - stack(i_stack - 1) = stack(i_stack - 1) .and. stack(i_stack) - i_stack = i_stack - 1 - case (OP_COMPLEMENT) - stack(i_stack) = .not. stack(i_stack) - case default - ! If the token is not an operator, evaluate the sense of particle with - ! respect to the surface and see if the token matches the sense. If the - ! particle's surface attribute is set and matches the token, that - ! overrides the determination based on sense(). - i_stack = i_stack + 1 - if (token == p%surface) then - stack(i_stack) = .true. - elseif (-token == p%surface) then - stack(i_stack) = .false. - else - actual_sense = surfaces(abs(token)) % sense( & - p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw) - stack(i_stack) = (actual_sense .eqv. (token > 0)) - end if - end select - - end do - - if (i_stack == 1) then - ! The one remaining logical on the stack indicates whether the particle is - ! in the cell. - in_cell = stack(i_stack) - else - ! This case occurs if there is no region specification since i_stack will - ! still be zero. - in_cell = .true. - end if - end function complex_cell_contains - !=============================================================================== ! CHECK_CELL_OVERLAP checks for overlapping cells at the current particle's ! position using cell_contains and the LocalCoord's built up by find_cell @@ -166,8 +77,8 @@ contains ! the particle should only be contained in one cell per level if (index_cell /= p % coord(j) % cell) then call fatal_error("Overlapping cells detected: " & - &// trim(to_str(cells(index_cell) % id)) // ", " & - &// trim(to_str(cells(p % coord(j) % cell) % id)) & + &// trim(to_str(cells(index_cell) % id())) // ", " & + &// trim(to_str(cells(p % coord(j) % cell) % id())) & &// " on universe " // trim(to_str(univ % id))) end if @@ -223,7 +134,7 @@ contains if (use_search_cells) then i_cell = search_cells(i) ! check to make sure search cell is in same universe - if (cells(i_cell) % universe /= i_universe) cycle + if (cells(i_cell) % universe() /= i_universe) cycle else i_cell = universes(i_universe) % cells(i) end if @@ -236,7 +147,7 @@ contains ! Show cell information on trace if (verbosity >= 10 .or. trace) then call write_message(" Entering cell " // trim(to_str(& - cells(i_cell) % id))) + cells(i_cell) % id()))) end if found = .true. @@ -246,7 +157,7 @@ contains if (found) then associate(c => cells(i_cell)) - CELL_TYPE: if (c % type == FILL_MATERIAL) then + CELL_TYPE: if (c % type() == FILL_MATERIAL) then ! ====================================================================== ! AT LOWEST UNIVERSE, TERMINATE SEARCH @@ -262,20 +173,20 @@ contains distribcell_index = c % distribcell_index offset = 0 do k = 1, p % n_coord - if (cells(p % coord(k) % cell) % type == FILL_UNIVERSE) then - offset = offset + cells(p % coord(k) % cell) % & - offset(distribcell_index) - elseif (cells(p % coord(k) % cell) % type == FILL_LATTICE) then + if (cells(p % coord(k) % cell) % type() == FILL_UNIVERSE) then + offset = offset + cells(p % coord(k) % cell) & + % offset(distribcell_index-1) + elseif (cells(p % coord(k) % cell) % type() == FILL_LATTICE) then if (lattices(p % coord(k + 1) % lattice) % obj & % are_valid_indices([& p % coord(k + 1) % lattice_x, & p % coord(k + 1) % lattice_y, & p % coord(k + 1) % lattice_z])) then - offset = offset + lattices(p % coord(k + 1) % lattice) % obj % & - offset(distribcell_index, & - p % coord(k + 1) % lattice_x, & + offset = offset + lattices(p % coord(k + 1) % lattice) % obj & + % offset(distribcell_index - 1, & + [p % coord(k + 1) % lattice_x, & p % coord(k + 1) % lattice_y, & - p % coord(k + 1) % lattice_z) + p % coord(k + 1) % lattice_z]) end if end if end do @@ -300,7 +211,7 @@ contains p % sqrtkT = c % sqrtkT(1) end if - elseif (c % type == FILL_UNIVERSE) then CELL_TYPE + elseif (c % type() == FILL_UNIVERSE) then CELL_TYPE ! ====================================================================== ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL @@ -311,7 +222,7 @@ contains ! Move particle to next level and set universe j = j + 1 p % n_coord = j - p % coord(j) % universe = c % fill + p % coord(j) % universe = c % fill() + 1 ! Apply translation if (allocated(c % translation)) then @@ -328,11 +239,11 @@ contains call find_cell(p, found) j = p % n_coord - elseif (c % type == FILL_LATTICE) then CELL_TYPE + elseif (c % type() == FILL_LATTICE) then CELL_TYPE ! ====================================================================== ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL - associate (lat => lattices(c % fill) % obj) + associate (lat => lattices(c % fill() + 1) % obj) ! Determine lattice indices i_xyz = lat % get_indices(p % coord(j) % xyz + TINY_BIT * p % coord(j) % uvw) @@ -341,7 +252,7 @@ contains p % coord(j + 1) % uvw = p % coord(j) % uvw ! set particle lattice indices - p % coord(j + 1) % lattice = c % fill + p % coord(j + 1) % lattice = c % fill() + 1 p % coord(j + 1) % lattice_x = i_xyz(1) p % coord(j + 1) % lattice_y = i_xyz(2) p % coord(j + 1) % lattice_z = i_xyz(3) @@ -350,18 +261,18 @@ contains if (lat % are_valid_indices(i_xyz)) then ! Particle is inside the lattice. p % coord(j + 1) % universe = & - lat % universes(i_xyz(1), i_xyz(2), i_xyz(3)) + lat % get([i_xyz(1), i_xyz(2), i_xyz(3)]) + 1 else ! Particle is outside the lattice. - if (lat % outer == NO_OUTER_UNIVERSE) then + if (lat % outer() == NO_OUTER_UNIVERSE) then call warning("Particle " // trim(to_str(p %id)) & - // " is outside lattice " // trim(to_str(lat % id)) & + // " is outside lattice " // trim(to_str(lat % id())) & // " but the lattice has no defined outer universe.") found = .false. return else - p % coord(j + 1) % universe = lat % outer + p % coord(j + 1) % universe = lat % outer() + 1 end if end if end associate @@ -396,7 +307,7 @@ contains lat => lattices(p % coord(j) % lattice) % obj if (verbosity >= 10 .or. trace) then - call write_message(" Crossing lattice " // trim(to_str(lat % id)) & + call write_message(" Crossing lattice " // trim(to_str(lat % id())) & &// ". Current position (" // trim(to_str(p % coord(j) % lattice_x)) & &// "," // trim(to_str(p % coord(j) % lattice_y)) // "," & &// trim(to_str(p % coord(j) % lattice_z)) // ")") @@ -428,7 +339,8 @@ contains else OUTSIDE_LAT ! Find cell in next lattice element - p % coord(j) % universe = lat % universes(i_xyz(1), i_xyz(2), i_xyz(3)) + p % coord(j) % universe = & + lat % get([i_xyz(1), i_xyz(2), i_xyz(3)]) + 1 call find_cell(p, found) if (.not. found) then @@ -465,26 +377,15 @@ contains integer, intent(out) :: lattice_translation(3) integer, intent(out) :: next_level - integer :: i ! index for surface in cell integer :: j - integer :: index_surf ! index in surfaces array (with sign) integer :: i_xyz(3) ! lattice indices integer :: level_surf_cross ! surface crossed on current level integer :: level_lat_trans(3) ! lattice translation on current level - real(8) :: x,y,z ! particle coordinates real(8) :: xyz_t(3) ! local particle coordinates - real(8) :: beta, gama ! skewed particle coordiantes - real(8) :: u,v,w ! particle directions - real(8) :: beta_dir ! skewed particle direction - real(8) :: gama_dir ! skewed particle direction - real(8) :: edge ! distance to oncoming edge - real(8) :: d ! evaluated distance real(8) :: d_lat ! distance to lattice boundary real(8) :: d_surf ! distance to surface - real(8) :: x0,y0,z0 ! coefficients for surface real(8) :: xyz_cross(3) ! coordinates at projected surface crossing real(8) :: surf_uvw(3) ! surface normal direction - logical :: coincident ! is particle on surface? type(Cell), pointer :: c class(Lattice), pointer :: lat @@ -502,35 +403,11 @@ contains ! get pointer to cell on this level c => cells(p % coord(j) % cell) - ! copy directional cosines - u = p % coord(j) % uvw(1) - v = p % coord(j) % uvw(2) - w = p % coord(j) % uvw(3) - ! ======================================================================= ! FIND MINIMUM DISTANCE TO SURFACE IN THIS CELL - SURFACE_LOOP: do i = 1, size(c % region) - index_surf = c % region(i) - coincident = (index_surf == p % surface) - - ! ignore this token if it corresponds to an operator rather than a - ! region. - index_surf = abs(index_surf) - if (index_surf >= OP_UNION) cycle - - ! Calculate distance to surface - d = surfaces(index_surf) % distance(p % coord(j) % xyz, & - p % coord(j) % uvw, logical(coincident, kind=C_BOOL)) - - ! Check if calculated distance is new minimum - if (d < d_surf) then - if (abs(d - d_surf)/d_surf >= FP_PRECISION) then - d_surf = d - level_surf_cross = -c % region(i) - end if - end if - end do SURFACE_LOOP + call c % distance(p % coord(j) % xyz, p % coord(j) % uvw, p % surface, & + d_surf, level_surf_cross) ! ======================================================================= ! FIND MINIMUM DISTANCE TO LATTICE SURFACES @@ -538,185 +415,22 @@ contains LAT_COORD: if (p % coord(j) % lattice /= NONE) then lat => lattices(p % coord(j) % lattice) % obj + i_xyz(1) = p % coord(j) % lattice_x + i_xyz(2) = p % coord(j) % lattice_y + i_xyz(3) = p % coord(j) % lattice_z + LAT_TYPE: select type(lat) type is (RectLattice) - ! copy local coordinates - x = p % coord(j) % xyz(1) - y = p % coord(j) % xyz(2) - z = p % coord(j) % xyz(3) - - ! determine oncoming edge - x0 = sign(lat % pitch(1) * HALF, u) - y0 = sign(lat % pitch(2) * HALF, v) - - ! left and right sides - if (abs(x - x0) < FP_PRECISION) then - d = INFINITY - elseif (u == ZERO) then - d = INFINITY - else - d = (x0 - x)/u - end if - - d_lat = d - if (u > 0) then - level_lat_trans(:) = [1, 0, 0] - else - level_lat_trans(:) = [-1, 0, 0] - end if - - ! front and back sides - if (abs(y - y0) < FP_PRECISION) then - d = INFINITY - elseif (v == ZERO) then - d = INFINITY - else - d = (y0 - y)/v - end if - - if (d < d_lat) then - d_lat = d - if (v > 0) then - level_lat_trans(:) = [0, 1, 0] - else - level_lat_trans(:) = [0, -1, 0] - end if - end if - - if (lat % is_3d) then - z0 = sign(lat % pitch(3) * HALF, w) - - ! top and bottom sides - if (abs(z - z0) < FP_PRECISION) then - d = INFINITY - elseif (w == ZERO) then - d = INFINITY - else - d = (z0 - z)/w - end if - - if (d < d_lat) then - d_lat = d - if (w > 0) then - level_lat_trans(:) = [0, 0, 1] - else - level_lat_trans(:) = [0, 0, -1] - end if - end if - end if + call lat % distance(p % coord(j) % xyz, p % coord(j) % uvw, & + i_xyz, d_lat, level_lat_trans) type is (HexLattice) LAT_TYPE - ! Copy local coordinates. - z = p % coord(j) % xyz(3) - i_xyz(1) = p % coord(j) % lattice_x - i_xyz(2) = p % coord(j) % lattice_y - i_xyz(3) = p % coord(j) % lattice_z - - ! Compute velocities along the hexagonal axes. - beta_dir = u*sqrt(THREE)/TWO + v/TWO - gama_dir = u*sqrt(THREE)/TWO - v/TWO - - ! Note that hexagonal lattice distance calculations are performed - ! using the particle's coordinates relative to the neighbor lattice - ! cells, not relative to the particle's current cell. This is done - ! because there is significant disagreement between neighboring cells - ! on where the lattice boundary is due to the worse finite precision - ! of hex lattices. - - ! Upper right and lower left sides. - edge = -sign(lat % pitch(1)/TWO, beta_dir) ! Oncoming edge - if (beta_dir > ZERO) then - xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[1, 0, 0]) - else - xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[-1, 0, 0]) - end if - beta = xyz_t(1)*sqrt(THREE)/TWO + xyz_t(2)/TWO - if (abs(beta - edge) < FP_PRECISION) then - d = INFINITY - else if (beta_dir == ZERO) then - d = INFINITY - else - d = (edge - beta)/beta_dir - end if - - d_lat = d - if (beta_dir > 0) then - level_lat_trans(:) = [1, 0, 0] - else - level_lat_trans(:) = [-1, 0, 0] - end if - - ! Lower right and upper left sides. - edge = -sign(lat % pitch(1)/TWO, gama_dir) ! Oncoming edge - if (gama_dir > ZERO) then - xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[1, -1, 0]) - else - xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[-1, 1, 0]) - end if - gama = xyz_t(1)*sqrt(THREE)/TWO - xyz_t(2)/TWO - if (abs(gama - edge) < FP_PRECISION) then - d = INFINITY - else if (gama_dir == ZERO) then - d = INFINITY - else - d = (edge - gama)/gama_dir - end if - - if (d < d_lat) then - d_lat = d - if (gama_dir > 0) then - level_lat_trans(:) = [1, -1, 0] - else - level_lat_trans(:) = [-1, 1, 0] - end if - end if - - ! Upper and lower sides. - edge = -sign(lat % pitch(1)/TWO, v) ! Oncoming edge - if (v > ZERO) then - xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[0, 1, 0]) - else - xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[0, -1, 0]) - end if - if (abs(xyz_t(2) - edge) < FP_PRECISION) then - d = INFINITY - else if (v == ZERO) then - d = INFINITY - else - d = (edge - xyz_t(2))/v - end if - - if (d < d_lat) then - d_lat = d - if (v > 0) then - level_lat_trans(:) = [0, 1, 0] - else - level_lat_trans(:) = [0, -1, 0] - end if - end if - - ! Top and bottom sides. - if (lat % is_3d) then - z0 = sign(lat % pitch(2) * HALF, w) - - if (abs(z - z0) < FP_PRECISION) then - d = INFINITY - elseif (w == ZERO) then - d = INFINITY - else - d = (z0 - z)/w - end if - - if (d < d_lat) then - d_lat = d - if (w > 0) then - level_lat_trans(:) = [0, 0, 1] - else - level_lat_trans(:) = [0, 0, -1] - end if - end if - end if + xyz_t(1) = p % coord(j-1) % xyz(1) + xyz_t(2) = p % coord(j-1) % xyz(2) + xyz_t(3) = p % coord(j) % xyz(3) + call lat % distance(xyz_t, p % coord(j) % uvw, & + i_xyz, d_lat, level_lat_trans) end select LAT_TYPE if (d_lat < ZERO) then @@ -738,7 +452,7 @@ contains ! positive half-space were given in the region specification. Thus, we ! have to explicitly check which half-space the particle would be ! traveling into if the surface is crossed - if (.not. c % simple) then + if (.not. c % simple()) then xyz_cross(:) = p % coord(j) % xyz + d_surf*p % coord(j) % uvw call surfaces(abs(level_surf_cross)) % normal(xyz_cross, surf_uvw) if (dot_product(p % coord(j) % uvw, surf_uvw) > ZERO) then @@ -820,388 +534,4 @@ contains end subroutine neighbor_lists -!=============================================================================== -! CALC_OFFSETS calculates and stores the offsets in all fill cells. This -! routine is called once upon initialization. -!=============================================================================== - - subroutine calc_offsets(univ_id, map, univ, counts, found) - - integer, intent(in) :: univ_id ! target universe ID - integer, intent(in) :: map ! map index in vector of maps - type(Universe), intent(in) :: univ ! universe searching in - integer, intent(inout) :: counts(:,:) ! target count - logical, intent(inout) :: found(:,:) ! target found - - integer :: i ! index over cells - integer :: j, k, m ! indices in lattice - integer :: n ! number of cells to search - integer :: offset ! total offset for a given cell - integer :: cell_index ! index in cells array - type(Cell), pointer :: c ! pointer to current cell - type(Universe), pointer :: next_univ ! next universe to cycle through - class(Lattice), pointer :: lat ! pointer to current lattice - - n = size(univ % cells) - offset = 0 - - do i = 1, n - - cell_index = univ % cells(i) - - ! get pointer to cell - c => cells(cell_index) - - ! ==================================================================== - ! AT LOWEST UNIVERSE, TERMINATE SEARCH - if (c % type == FILL_MATERIAL) then - - ! ==================================================================== - ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL - elseif (c % type == FILL_UNIVERSE) then - ! Set offset for the cell on this level - c % offset(map) = offset - - ! Count contents of this cell - next_univ => universes(c % fill) - offset = offset + count_target(next_univ, counts, found, univ_id, map) - - ! Move into the next universe - next_univ => universes(c % fill) - c => cells(cell_index) - - ! ==================================================================== - ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL - elseif (c % type == FILL_LATTICE) then - - ! Set current lattice - lat => lattices(c % fill) % obj - - select type (lat) - - type is (RectLattice) - - ! Loop over lattice coordinates - do j = 1, lat % n_cells(1) - do k = 1, lat % n_cells(2) - do m = 1, lat % n_cells(3) - lat % offset(map, j, k, m) = offset - next_univ => universes(lat % universes(j, k, m)) - offset = offset + & - count_target(next_univ, counts, found, univ_id, map) - end do - end do - end do - - type is (HexLattice) - - ! Loop over lattice coordinates - do m = 1, lat % n_axial - do k = 1, 2*lat % n_rings - 1 - do j = 1, 2*lat % n_rings - 1 - ! This array location is never used - if (j + k < lat % n_rings + 1) then - cycle - ! This array location is never used - else if (j + k > 3*lat % n_rings - 1) then - cycle - else - lat % offset(map, j, k, m) = offset - next_univ => universes(lat % universes(j, k, m)) - offset = offset + & - count_target(next_univ, counts, found, univ_id, map) - end if - end do - end do - end do - end select - - end if - end do - - end subroutine calc_offsets - -!=============================================================================== -! COUNT_TARGET recursively totals the numbers of occurances of a given -! universe ID beginning with the universe given. -!=============================================================================== - - recursive function count_target(univ, counts, found, univ_id, map) result(count) - - type(Universe), intent(in) :: univ ! universe to search through - integer, intent(inout) :: counts(:,:) ! target count - logical, intent(inout) :: found(:,:) ! target found - integer, intent(in) :: univ_id ! target universe ID - integer, intent(in) :: map ! current map - - integer :: i ! index over cells - integer :: j, k, m ! indices in lattice - integer :: n ! number of cells to search - integer :: cell_index ! index in cells array - integer :: count ! number of times target located - type(Cell), pointer :: c ! pointer to current cell - type(Universe), pointer :: next_univ ! next univ to loop through - class(Lattice), pointer :: lat ! pointer to current lattice - - ! Don't research places already checked - if (found(universe_dict % get(univ % id), map)) then - count = counts(universe_dict % get(univ % id), map) - return - end if - - ! If this is the target, it can't contain itself. - ! Count = 1, then quit - if (univ % id == univ_id) then - count = 1 - counts(universe_dict % get(univ % id), map) = 1 - found(universe_dict % get(univ % id), map) = .true. - return - end if - - count = 0 - n = size(univ % cells) - - do i = 1, n - - cell_index = univ % cells(i) - - ! get pointer to cell - c => cells(cell_index) - - ! ==================================================================== - ! AT LOWEST UNIVERSE, TERMINATE SEARCH - if (c % type == FILL_MATERIAL) then - - ! ==================================================================== - ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL - elseif (c % type == FILL_UNIVERSE) then - - next_univ => universes(c % fill) - - ! Found target - stop since target cannot contain itself - if (next_univ % id == univ_id) then - count = count + 1 - return - end if - - count = count + count_target(next_univ, counts, found, univ_id, map) - c => cells(cell_index) - - ! ==================================================================== - ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL - elseif (c % type == FILL_LATTICE) then - - ! Set current lattice - lat => lattices(c % fill) % obj - - select type (lat) - - type is (RectLattice) - - ! Loop over lattice coordinates - do j = 1, lat % n_cells(1) - do k = 1, lat % n_cells(2) - do m = 1, lat % n_cells(3) - next_univ => universes(lat % universes(j, k, m)) - - ! Found target - stop since target cannot contain itself - if (next_univ % id == univ_id) then - count = count + 1 - cycle - end if - - count = count + & - count_target(next_univ, counts, found, univ_id, map) - - end do - end do - end do - - type is (HexLattice) - - ! Loop over lattice coordinates - do m = 1, lat % n_axial - do k = 1, 2*lat % n_rings - 1 - do j = 1, 2*lat % n_rings - 1 - ! This array location is never used - if (j + k < lat % n_rings + 1) then - cycle - ! This array location is never used - else if (j + k > 3*lat % n_rings - 1) then - cycle - else - next_univ => universes(lat % universes(j, k, m)) - - ! Found target - stop since target cannot contain itself - if (next_univ % id == univ_id) then - count = count + 1 - cycle - end if - - count = count + & - count_target(next_univ, counts, found, univ_id, map) - end if - end do - end do - end do - - end select - - end if - end do - - counts(universe_dict % get(univ % id), map) = count - found(universe_dict % get(univ % id), map) = .true. - - end function count_target - -!=============================================================================== -! COUNT_INSTANCE recursively totals the number of occurrences of all cells -! beginning with the universe given. -!=============================================================================== - - recursive subroutine count_instance(univ) - - type(Universe), intent(in) :: univ ! universe to search through - - integer :: i ! index over cells - integer :: j, k, m ! indices in lattice - integer :: n ! number of cells to search - - n = size(univ % cells) - - do i = 1, n - associate (c => cells(univ % cells(i))) - c % instances = c % instances + 1 - - ! ==================================================================== - ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL - if (c % type == FILL_UNIVERSE) then - - call count_instance(universes(c % fill)) - - ! ==================================================================== - ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL - elseif (c % type == FILL_LATTICE) then - - ! Set current lattice - associate (lat => lattices(c % fill) % obj) - - select type (lat) - type is (RectLattice) - - ! Loop over lattice coordinates - do j = 1, lat % n_cells(1) - do k = 1, lat % n_cells(2) - do m = 1, lat % n_cells(3) - call count_instance(universes(lat % universes(j, k, m))) - end do - end do - end do - - type is (HexLattice) - - ! Loop over lattice coordinates - do m = 1, lat % n_axial - do k = 1, 2*lat % n_rings - 1 - do j = 1, 2*lat % n_rings - 1 - ! This array location is never used - if (j + k < lat % n_rings + 1) then - cycle - ! This array location is never used - else if (j + k > 3*lat % n_rings - 1) then - cycle - else - call count_instance(universes(lat % universes(j, k, m))) - end if - end do - end do - end do - - end select - end associate - end if - end associate - end do - - end subroutine count_instance - -!=============================================================================== -! MAXIMUM_LEVELS determines the maximum number of nested coordinate levels in -! the geometry -!=============================================================================== - - recursive function maximum_levels(univ) result(levels) - - type(Universe), intent(in) :: univ ! universe to search through - integer :: levels ! maximum number of levels for this universe - - integer :: i ! index over cells - integer :: j, k, m ! indices in lattice - integer :: levels_below ! max levels below this universe - type(Cell), pointer :: c ! pointer to current cell - type(Universe), pointer :: next_univ ! next universe to loop through - class(Lattice), pointer :: lat ! pointer to current lattice - - levels_below = 0 - do i = 1, size(univ % cells) - c => cells(univ % cells(i)) - - ! ==================================================================== - ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL - if (c % type == FILL_UNIVERSE) then - - next_univ => universes(c % fill) - levels_below = max(levels_below, maximum_levels(next_univ)) - - ! ==================================================================== - ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL - elseif (c % type == FILL_LATTICE) then - - ! Set current lattice - lat => lattices(c % fill) % obj - - select type (lat) - - type is (RectLattice) - - ! Loop over lattice coordinates - do j = 1, lat % n_cells(1) - do k = 1, lat % n_cells(2) - do m = 1, lat % n_cells(3) - next_univ => universes(lat % universes(j, k, m)) - levels_below = max(levels_below, maximum_levels(next_univ)) - end do - end do - end do - - type is (HexLattice) - - ! Loop over lattice coordinates - do m = 1, lat % n_axial - do k = 1, 2*lat % n_rings - 1 - do j = 1, 2*lat % n_rings - 1 - ! This array location is never used - if (j + k < lat % n_rings + 1) then - cycle - ! This array location is never used - else if (j + k > 3*lat % n_rings - 1) then - cycle - else - next_univ => universes(lat % universes(j, k, m)) - levels_below = max(levels_below, maximum_levels(next_univ)) - end if - end do - end do - end do - - end select - - end if - end do - - levels = 1 + levels_below - - end function maximum_levels - end module geometry diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp new file mode 100644 index 0000000000..3f36cf299a --- /dev/null +++ b/src/geometry_aux.cpp @@ -0,0 +1,341 @@ +#include "geometry_aux.h" + +#include // for std::max +#include +#include + +#include "cell.h" +#include "constants.h" +#include "error.h" +#include "lattice.h" + + +namespace openmc { + +//============================================================================== + +void +adjust_indices_c() +{ + // Adjust material/fill idices. + for (Cell *c : global_cells) { + if (c->material[0] == C_NONE) { + int32_t id = c->fill; + auto search_univ = universe_map.find(id); + auto search_lat = lattice_map.find(id); + if (search_univ != universe_map.end()) { + c->type = FILL_UNIVERSE; + c->fill = search_univ->second; + } else if (search_lat != lattice_map.end()) { + c->type = FILL_LATTICE; + c->fill = search_lat->second; + } else { + std::stringstream err_msg; + err_msg << "Specified fill " << id << " on cell " << c->id + << " is neither a universe nor a lattice."; + fatal_error(err_msg); + } + } else { + //TODO: materials + c->type = FILL_MATERIAL; + } + } + + // Change cell.universe values from IDs to indices. + for (Cell *c : global_cells) { + auto search = universe_map.find(c->universe); + if (search != universe_map.end()) { + //TODO: Remove this off-by-one indexing. + c->universe = search->second + 1; + } else { + std::stringstream err_msg; + err_msg << "Could not find universe " << c->universe + << " specified on cell " << c->id; + fatal_error(err_msg); + } + } + + // Change all lattice universe values from IDs to indices. + for (Lattice *l : lattices_c) { + l->adjust_indices(); + } +} + +//============================================================================== + +int32_t +find_root_universe() +{ + // Find all the universes listed as a cell fill. + std::unordered_set fill_univ_ids; + for (Cell *c : global_cells) { + fill_univ_ids.insert(c->fill); + } + + // Find all the universes contained in a lattice. + for (Lattice *lat : lattices_c) { + for (auto it = lat->begin(); it != lat->end(); ++it) { + fill_univ_ids.insert(*it); + } + if (lat->outer != NO_OUTER_UNIVERSE) { + fill_univ_ids.insert(lat->outer); + } + } + + // Figure out which universe is not in the set. This is the root universe. + bool root_found {false}; + int32_t root_univ; + for (int32_t i = 0; i < global_universes.size(); i++) { + auto search = fill_univ_ids.find(global_universes[i]->id); + if (search == fill_univ_ids.end()) { + if (root_found) { + fatal_error("Two or more universes are not used as fill universes, so " + "it is not possible to distinguish which one is the root " + "universe."); + } else { + root_found = true; + root_univ = i; + } + } + } + if (!root_found) fatal_error("Could not find a root universe. Make sure " + "there are no circular dependencies in the geometry."); + + return root_univ; +} + +//============================================================================== + +void +allocate_offset_tables(int n_maps) +{ + for (Cell *c : global_cells) { + if (c->type != FILL_MATERIAL) { + c->offset.resize(n_maps, C_NONE); + } + } + + for (Lattice *lat : lattices_c) { + lat->allocate_offset_table(n_maps); + } +} + +//============================================================================== + +void +count_cell_instances(int32_t univ_indx) +{ + for (int32_t cell_indx : global_universes[univ_indx]->cells) { + Cell &c = *global_cells[cell_indx]; + ++c.n_instances; + + if (c.type == FILL_UNIVERSE) { + // This cell contains another universe. Recurse into that universe. + count_cell_instances(c.fill); + + } else if (c.type == FILL_LATTICE) { + // This cell contains a lattice. Recurse into the lattice universes. + Lattice &lat = *lattices_c[c.fill]; + for (auto it = lat.begin(); it != lat.end(); ++it) { + count_cell_instances(*it); + } + } + } +} + +//============================================================================== + +int +count_universe_instances(int32_t search_univ, int32_t target_univ_id) +{ + // If this is the target, it can't contain itself. + if (global_universes[search_univ]->id == target_univ_id) { + return 1; + } + + int count {0}; + for (int32_t cell_indx : global_universes[search_univ]->cells) { + Cell &c = *global_cells[cell_indx]; + + if (c.type == FILL_UNIVERSE) { + int32_t next_univ = c.fill; + count += count_universe_instances(next_univ, target_univ_id); + + } else if (c.type == FILL_LATTICE) { + Lattice &lat = *lattices_c[c.fill]; + for (auto it = lat.begin(); it != lat.end(); ++it) { + int32_t next_univ = *it; + count += count_universe_instances(next_univ, target_univ_id); + } + } + } + + return count; +} + +//============================================================================== + +void +fill_offset_tables(int32_t target_univ_id, int map) +{ + for (Universe *univ : global_universes) { + int32_t offset {0}; // TODO: is this a bug? It matches F90 implementation. + for (int32_t cell_indx : univ->cells) { + Cell &c = *global_cells[cell_indx]; + + if (c.type == FILL_UNIVERSE) { + c.offset[map] = offset; + int32_t search_univ = c.fill; + offset += count_universe_instances(search_univ, target_univ_id); + + } else if (c.type == FILL_LATTICE) { + Lattice &lat = *lattices_c[c.fill]; + offset = lat.fill_offset_table(offset, target_univ_id, map); + } + } + } +} + +//============================================================================== + +std::string +distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset, + const Universe &search_univ, int32_t offset) +{ + std::stringstream path; + + path << "u" << search_univ.id << "->"; + + // Check to see if this universe directly contains the target cell. If so, + // write to the path and return. + for (int32_t cell_indx : search_univ.cells) { + if ((cell_indx == target_cell) && (offset == target_offset)) { + Cell &c = *global_cells[cell_indx]; + path << "c" << c.id; + return path.str(); + } + } + + // The target must be further down the geometry tree and contained in a fill + // cell or lattice cell in this universe. Find which cell contains the + // target. + std::vector::const_reverse_iterator cell_it + {search_univ.cells.crbegin()}; + for (; cell_it != search_univ.cells.crend(); ++cell_it) { + Cell &c = *global_cells[*cell_it]; + + // Material cells don't contain other cells so ignore them. + if (c.type != FILL_MATERIAL) { + int32_t temp_offset; + if (c.type == FILL_UNIVERSE) { + temp_offset = offset + c.offset[map]; + } else { + Lattice &lat = *lattices_c[c.fill]; + int32_t indx = lat.universes.size()*map + lat.begin().indx; + temp_offset = offset + lat.offsets[indx]; + } + + // The desired cell is the first cell that gives an offset smaller or + // equal to the target offset. + if (temp_offset <= target_offset) break; + } + } + + // Add the cell to the path string. + Cell &c = *global_cells[*cell_it]; + path << "c" << c.id << "->"; + + if (c.type == FILL_UNIVERSE) { + // Recurse into the fill cell. + offset += c.offset[map]; + path << distribcell_path_inner(target_cell, map, target_offset, + *global_universes[c.fill], offset); + return path.str(); + } else { + // Recurse into the lattice cell. + Lattice &lat = *lattices_c[c.fill]; + path << "l" << lat.id; + for (ReverseLatticeIter it = lat.rbegin(); it != lat.rend(); ++it) { + int32_t indx = lat.universes.size()*map + it.indx; + int32_t temp_offset = offset + lat.offsets[indx]; + if (temp_offset <= target_offset) { + offset = temp_offset; + path << "(" << lat.index_to_string(it.indx) << ")->"; + path << distribcell_path_inner(target_cell, map, target_offset, + *global_universes[*it], offset); + return path.str(); + } + } + } +} + +//============================================================================== + +int +distribcell_path_len(int32_t target_cell, int32_t map, int32_t target_offset, + int32_t root_univ) +{ + Universe &root = *global_universes[root_univ]; + std::string path_ {distribcell_path_inner(target_cell, map, target_offset, + root, 0)}; + return path_.size() + 1; +} + +//============================================================================== + +void +distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset, + int32_t root_univ, char *path) +{ + Universe &root = *global_universes[root_univ]; + std::string path_ {distribcell_path_inner(target_cell, map, target_offset, + root, 0)}; + path_.copy(path, path_.size()); + path[path_.size()] = '\0'; +} + +//============================================================================== + +int +maximum_levels(int32_t univ) +{ + int levels_below {0}; + + for (int32_t cell_indx : global_universes[univ]->cells) { + Cell &c = *global_cells[cell_indx]; + if (c.type == FILL_UNIVERSE) { + int32_t next_univ = c.fill; + levels_below = std::max(levels_below, maximum_levels(next_univ)); + } else if (c.type == FILL_LATTICE) { + Lattice &lat = *lattices_c[c.fill]; + for (auto it = lat.begin(); it != lat.end(); ++it) { + int32_t next_univ = *it; + levels_below = std::max(levels_below, maximum_levels(next_univ)); + } + } + } + + ++levels_below; + return levels_below; +} + +//============================================================================== + +void +free_memory_geometry_c() +{ + for (Cell *c : global_cells) {delete c;} + global_cells.clear(); + cell_map.clear(); + n_cells = 0; + + for (Universe *u : global_universes) {delete u;} + global_universes.clear(); + universe_map.clear(); + + for (Lattice *lat : lattices_c) {delete lat;} + lattices_c.clear(); + lattice_map.clear(); +} + +} // namespace openmc diff --git a/src/geometry_aux.h b/src/geometry_aux.h new file mode 100644 index 0000000000..b998ad32c9 --- /dev/null +++ b/src/geometry_aux.h @@ -0,0 +1,112 @@ +//! \file geometry_aux.h +//! Auxilary functions for geometry initialization and general data handling. + +#ifndef GEOMETRY_AUX_H +#define GEOMETRY_AUX_H + +#include + + +namespace openmc { + +//============================================================================== +//! Replace Universe, Lattice, and Material IDs with indices. +//============================================================================== + +extern "C" void adjust_indices_c(); + +//============================================================================== +//! Figure out which Universe is the root universe. +//! +//! This function looks for a universe that is not listed in a Cell::fill or in +//! a Lattice. +//! @return The index of the root universe. +//============================================================================== + +extern "C" int32_t find_root_universe(); + +//============================================================================== +//! Allocate storage in Lattice and Cell objects for distribcell offset tables. +//============================================================================== + +extern "C" void allocate_offset_tables(int n_maps); + +//============================================================================== +//! Recursively search through the geometry and count cell instances. +//! +//! This function will update the Cell::n_instances value for each cell in the +//! geometry. +//! @param univ_indx The index of the universe to begin searching from (probably +//! the root universe). +//============================================================================== + +extern "C" void count_cell_instances(int32_t univ_indx); + +//============================================================================== +//! Recursively search through universes and count universe instances. +//! @param search_univ The index of the universe to begin searching from. +//! @param target_univ_id The ID of the universe to be counted. +//! @return The number of instances of target_univ_id in the geometry tree under +//! search_univ. +//============================================================================== + +extern "C" int +count_universe_instances(int32_t search_univ, int32_t target_univ_id); + +//============================================================================== +//! Populate Cell and Lattice distribcell offset tables. +//! @param target_univ_id The ID of the universe to be counted. +//! @param map The index of the distribcell map that defines the offsets for the +//! target universe. +//============================================================================== + +extern "C" void fill_offset_tables(int32_t target_univ_id, int map); + +//============================================================================== +//! Find the length necessary for a string to contain a distribcell path. +//! @param target_cell The index of the Cell in the global Cell array. +//! @param map The index of the distribcell mapping corresponding to the target +//! cell. +//! @param target_offset An instance number for a distributed cell. +//! @param root_univ The index of the root Universe in the global Universe +//! array. +//! @return The size of a character array needed to fit the distribcell path. +//============================================================================== + +extern "C" int +distribcell_path_len(int32_t target_cell, int32_t map, int32_t target_offset, + int32_t root_univ); + +//============================================================================== +//! Build a character array representing the path to a distribcell instance. +//! @param target_cell The index of the Cell in the global Cell array. +//! @param map The index of the distribcell mapping corresponding to the target +//! cell. +//! @param target_offset An instance number for a distributed cell. +//! @param root_univ The index of the root Universe in the global Universe +//! array. +//! @param[out] path The unique traversal through the geometry tree that leads +//! to the desired instance of the target cell. +//============================================================================== + +extern "C" void +distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset, + int32_t root_univ, char *path); + +//============================================================================== +//! Determine the maximum number of nested coordinate levels in the geometry. +//! @param univ The index of the universe to begin seraching from (probably the +//! root universe). +//! @return The number of coordinate levels. +//============================================================================== + +extern "C" int maximum_levels(int32_t univ); + +//============================================================================== +//! Deallocates global vectors and maps for cells, universes, and lattices. +//============================================================================== + +extern "C" void free_memory_geometry_c(); + +} // namespace openmc +#endif // GEOMETRY_AUX_H diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index f6206adf99..758fc54da7 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -6,6 +6,7 @@ module geometry_header use constants, only: HALF, TWO, THREE, INFINITY, K_BOLTZMANN, & MATERIAL_VOID, NONE use dict_header, only: DictCharInt, DictIntInt + use hdf5_interface, only: HID_T use material_header, only: Material, materials, material_dict, n_materials use nuclide_header use sab_header @@ -14,13 +15,192 @@ module geometry_header implicit none + interface + function cell_pointer_c(cell_ind) bind(C, name='cell_pointer') result(ptr) + import C_PTR, C_INT32_T + integer(C_INT32_T), intent(in), value :: cell_ind + type(C_PTR) :: ptr + end function cell_pointer_c + + function cell_id_c(cell_ptr) bind(C, name='cell_id') result(id) + import C_PTR, C_INT32_T + type(C_PTR), intent(in), value :: cell_ptr + integer(C_INT32_T) :: id + end function cell_id_c + + subroutine cell_set_id_c(cell_ptr, id) bind(C, name='cell_set_id') + import C_PTR, C_INT32_T + type(C_PTR), intent(in), value :: cell_ptr + integer(C_INT32_T), intent(in), value :: id + end subroutine cell_set_id_c + + function cell_type_c(cell_ptr) bind(C, name='cell_type') result(type) + import C_PTR, C_INT + type(C_PTR), intent(in), value :: cell_ptr + integer(C_INT) :: type + end function cell_type_c + + subroutine cell_set_type_c(cell_ptr, type) bind(C, name='cell_set_type') + import C_PTR, C_INT + type(C_PTR), intent(in), value :: cell_ptr + integer(C_INT), intent(in), value :: type + end subroutine cell_set_type_c + + function cell_universe_c(cell_ptr) bind(C, name='cell_universe') & + result(universe) + import C_PTR, C_INT32_T + type(C_PTR), intent(in), value :: cell_ptr + integer(C_INT32_T) :: universe + end function cell_universe_c + + subroutine cell_set_universe_c(cell_ptr, universe) & + bind(C, name='cell_set_universe') + import C_PTR, C_INT32_T + type(C_PTR), intent(in), value :: cell_ptr + integer(C_INT32_T), intent(in), value :: universe + end subroutine cell_set_universe_c + + function cell_fill_c(cell_ptr) bind(C, name="cell_fill") result(fill) + import C_PTR, C_INT32_T + type(C_PTR), intent(in), value :: cell_ptr + integer(C_INT32_T) :: fill + end function cell_fill_c + + function cell_fill_ptr(cell_ptr) bind(C) result(fill_ptr) + import C_PTR + type(C_PTR), intent(in), value :: cell_ptr + type(C_PTR) :: fill_ptr + end function cell_fill_ptr + + function cell_n_instances_c(cell_ptr) bind(C, name='cell_n_instances') & + result(n_instances) + import C_PTR, C_INT32_T + type(C_PTR), intent(in), value :: cell_ptr + integer(C_INT32_T) :: n_instances + end function cell_n_instances_c + + function cell_simple_c(cell_ptr) bind(C, name='cell_simple') result(simple) + import C_PTR, C_BOOL + type(C_PTR), intent(in), value :: cell_ptr + logical(C_BOOL) :: simple + end function cell_simple_c + + subroutine cell_distance_c(cell_ptr, xyz, uvw, on_surface, min_dist, & + i_surf) bind(C, name="cell_distance") + import C_PTR, C_INT32_T, C_DOUBLE + type(C_PTR), intent(in), value :: cell_ptr + real(C_DOUBLE), intent(in) :: xyz(3) + real(C_DOUBLE), intent(in) :: uvw(3) + integer(C_INT32_T), intent(in), value :: on_surface + real(C_DOUBLE), intent(out) :: min_dist + integer(C_INT32_T), intent(out) :: i_surf + end subroutine cell_distance_c + + function cell_offset_c(cell_ptr, map) bind(C, name="cell_offset") & + result(offset) + import C_PTR, C_INT, C_INT32_T + type(C_PTR), intent(in), value :: cell_ptr + integer(C_INT), intent(in), value :: map + integer(C_INT32_T) :: offset + end function cell_offset_c + + subroutine cell_to_hdf5_c(cell_ptr, group) bind(C, name='cell_to_hdf5') + import HID_T, C_PTR + type(C_PTR), intent(in), value :: cell_ptr + integer(HID_T), intent(in), value :: group + end subroutine cell_to_hdf5_c + + function lattice_pointer_c(lat_ind) bind(C, name='lattice_pointer') & + result(ptr) + import C_PTR, C_INT32_T + integer(C_INT32_T), intent(in), value :: lat_ind + type(C_PTR) :: ptr + end function lattice_pointer_c + + function lattice_id_c(lat_ptr) bind(C, name='lattice_id') result(id) + import C_PTR, C_INT32_T + type(C_PTR), intent(in), value :: lat_ptr + integer(C_INT32_T) :: id + end function lattice_id_c + + function lattice_are_valid_indices_c(lat_ptr, i_xyz) & + bind(C, name='lattice_are_valid_indices') result (is_valid) + import C_PTR, C_INT, C_BOOL + type(C_PTR), intent(in), value :: lat_ptr + integer(C_INT), intent(in) :: i_xyz(3) + logical(C_BOOL) :: is_valid + end function lattice_are_valid_indices_c + + subroutine lattice_distance_c(lat_ptr, xyz, uvw, i_xyz, d, lattice_trans) & + bind(C, name='lattice_distance') + import C_PTR, C_INT, C_DOUBLE + type(C_PTR), intent(in), value :: lat_ptr + real(C_DOUBLE), intent(in) :: xyz(3) + real(C_DOUBLE), intent(in) :: uvw(3) + integer(C_INT), intent(in) :: i_xyz(3) + real(C_DOUBLE), intent(out) :: d + integer(C_INT), intent(out) :: lattice_trans(3) + end subroutine lattice_distance_c + + subroutine lattice_get_indices_c(lat_ptr, xyz, i_xyz) & + bind(C, name='lattice_get_indices') + import C_PTR, C_INT, C_DOUBLE + type(C_PTR), intent(in), value :: lat_ptr + real(C_DOUBLE), intent(in) :: xyz(3) + integer(C_INT), intent(out) :: i_xyz(3) + end subroutine lattice_get_indices_c + + subroutine lattice_get_local_xyz_c(lat_ptr, global_xyz, i_xyz, local_xyz) & + bind(C, name='lattice_get_local_xyz') + import C_PTR, C_INT, C_DOUBLE + type(C_PTR), intent(in), value :: lat_ptr + real(C_DOUBLE), intent(in) :: global_xyz(3) + integer(C_INT), intent(in) :: i_xyz(3) + real(C_DOUBLE), intent(out) :: local_xyz(3) + end subroutine lattice_get_local_xyz_c + + subroutine lattice_to_hdf5_c(lat_ptr, group) bind(C, name='lattice_to_hdf5') + import HID_T, C_PTR + type(C_PTR), intent(in), value :: lat_ptr + integer(HID_T), intent(in), value :: group + end subroutine lattice_to_hdf5_c + + function lattice_offset_c(lat_ptr, map, i_xyz) & + bind(C, name='lattice_offset') result(offset) + import C_PTR, C_INT, C_INT32_T + type(C_PTR), intent(in), value :: lat_ptr + integer(C_INT), intent(in), value :: map + integer(C_INT), intent(in) :: i_xyz(3) + integer(C_INT32_T) :: offset + end function lattice_offset_c + + function lattice_outer_c(lat_ptr) bind(C, name='lattice_outer') & + result(outer) + import C_PTR, C_INT32_T + type(C_PTR), intent(in), value :: lat_ptr + integer(C_INT32_T) :: outer + end function lattice_outer_c + + function lattice_universe_c(lat_ptr, i_xyz) & + bind(C, name='lattice_universe') result(univ) + import C_PTR, C_INT32_t, C_INT + type(C_PTR), intent(in), value :: lat_ptr + integer(C_INT), intent(in) :: i_xyz(3) + integer(C_INT32_T) :: univ + end function lattice_universe_c + + subroutine extend_cells_c(n) bind(C) + import C_INT32_t + integer(C_INT32_T), intent(in), value :: n + end subroutine extend_cells_c + end interface + !=============================================================================== ! UNIVERSE defines a geometry that fills all phase space !=============================================================================== type Universe integer :: id ! Unique ID - integer :: type ! Type integer, allocatable :: cells(:) ! List of cells within real(8) :: x0 ! Translation in x-coordinate real(8) :: y0 ! Translation in y-coordinate @@ -32,68 +212,24 @@ module geometry_header !=============================================================================== type, abstract :: Lattice - integer :: id ! Universe number for lattice - character(len=104) :: name = "" ! User-defined name - real(8), allocatable :: pitch(:) ! Pitch along each axis - integer, allocatable :: universes(:,:,:) ! Specified universes - integer :: outside ! Material to fill area outside - integer :: outer ! universe to tile outside the lat - logical :: is_3d ! Lattice has cells on z axis - integer, allocatable :: offset(:,:,:,:) ! Distribcell offsets + type(C_PTR) :: ptr contains - procedure(lattice_are_valid_indices_), deferred :: are_valid_indices - procedure(lattice_get_indices_), deferred :: get_indices - procedure(lattice_get_local_xyz_), deferred :: get_local_xyz + procedure :: id => lattice_id + procedure :: are_valid_indices => lattice_are_valid_indices + procedure :: distance => lattice_distance + procedure :: get => lattice_get + procedure :: get_indices => lattice_get_indices + procedure :: get_local_xyz => lattice_get_local_xyz + procedure :: offset => lattice_offset + procedure :: outer => lattice_outer + procedure :: to_hdf5 => lattice_to_hdf5 end type Lattice - abstract interface - -!=============================================================================== -! ARE_VALID_INDICES returns .true. if the given lattice indices fit within the -! bounds of the lattice. Returns false otherwise. - - function lattice_are_valid_indices_(this, i_xyz) result(is_valid) - import Lattice - class(Lattice), intent(in) :: this - integer, intent(in) :: i_xyz(3) - logical :: is_valid - end function lattice_are_valid_indices_ - -!=============================================================================== -! GET_INDICES returns the indices in a lattice for the given global xyz. - - function lattice_get_indices_(this, global_xyz) result(i_xyz) - import Lattice - class(Lattice), intent(in) :: this - real(8), intent(in) :: global_xyz(3) - integer :: i_xyz(3) - end function lattice_get_indices_ - -!=============================================================================== -! GET_LOCAL_XYZ returns the translated local version of the given global xyz. - - function lattice_get_local_xyz_(this, global_xyz, i_xyz) result(local_xyz) - import Lattice - class(Lattice), intent(in) :: this - real(8), intent(in) :: global_xyz(3) - integer, intent(in) :: i_xyz(3) - real(8) :: local_xyz(3) - end function lattice_get_local_xyz_ - end interface - !=============================================================================== ! RECTLATTICE extends LATTICE for rectilinear arrays. !=============================================================================== type, extends(Lattice) :: RectLattice - integer :: n_cells(3) ! Number of cells along each axis - real(8), allocatable :: lower_left(:) ! Global lower-left corner of lat - - contains - - procedure :: are_valid_indices => valid_inds_rect - procedure :: get_indices => get_inds_rect - procedure :: get_local_xyz => get_local_rect end type RectLattice !=============================================================================== @@ -101,15 +237,6 @@ module geometry_header !=============================================================================== type, extends(Lattice) :: HexLattice - integer :: n_rings ! Number of radial ring cell positoins - integer :: n_axial ! Number of axial cell positions - real(8), allocatable :: center(:) ! Global center of lattice - - contains - - procedure :: are_valid_indices => valid_inds_hex - procedure :: get_indices => get_inds_hex - procedure :: get_local_xyz => get_local_hex end type HexLattice !=============================================================================== @@ -125,25 +252,13 @@ module geometry_header !=============================================================================== type Cell - integer :: id ! Unique ID - character(len=104) :: name = "" ! User-defined name - integer :: type ! Type of cell (normal, universe, - ! lattice) - integer :: universe ! universe # this cell is in - integer :: fill ! universe # filling this cell - integer :: instances ! number of instances of this cell in - ! the geom + type(C_PTR) :: ptr + integer, allocatable :: material(:) ! Material within cell. Multiple ! materials for distribcell ! instances. 0 signifies a universe - integer, allocatable :: offset(:) ! Distribcell offset for tally - ! counter integer, allocatable :: region(:) ! Definition of spatial region as ! Boolean expression of half-spaces - integer, allocatable :: rpn(:) ! Reverse Polish notation for region - ! expression - logical :: simple ! Is the region simple (intersections - ! only) integer :: distribcell_index ! Index corresponding to this cell in ! distribcell arrays real(8), allocatable :: sqrtkT(:) ! Square root of k_Boltzmann * @@ -154,6 +269,22 @@ module geometry_header real(8), allocatable :: translation(:) real(8), allocatable :: rotation(:) real(8), allocatable :: rotation_matrix(:,:) + + contains + + procedure :: id => cell_id + procedure :: set_id => cell_set_id + procedure :: type => cell_type + procedure :: set_type => cell_set_type + procedure :: universe => cell_universe + procedure :: set_universe => cell_set_universe + procedure :: fill => cell_fill + procedure :: n_instances => cell_n_instances + procedure :: simple => cell_simple + procedure :: distance => cell_distance + procedure :: offset => cell_offset + procedure :: to_hdf5 => cell_to_hdf5 + end type Cell ! array index of the root universe @@ -161,7 +292,6 @@ module geometry_header integer(C_INT32_T), bind(C) :: n_cells ! # of cells integer(C_INT32_T), bind(C) :: n_universes ! # of universes - integer(C_INT32_T), bind(C) :: n_lattices ! # of lattices type(Cell), allocatable, target :: cells(:) type(Universe), allocatable, target :: universes(:) @@ -174,172 +304,150 @@ module geometry_header contains -!=============================================================================== + function lattice_id(this) result(id) + class(Lattice), intent(in) :: this + integer(C_INT32_T) :: id + id = lattice_id_c(this % ptr) + end function lattice_id - function valid_inds_rect(this, i_xyz) result(is_valid) - class(RectLattice), intent(in) :: this - integer, intent(in) :: i_xyz(3) - logical :: is_valid + function lattice_are_valid_indices(this, i_xyz) result (is_valid) + class(Lattice), intent(in) :: this + integer(C_INT), intent(in) :: i_xyz(3) + logical(C_BOOL) :: is_valid + is_valid = lattice_are_valid_indices_c(this % ptr, i_xyz) + end function lattice_are_valid_indices - is_valid = all(i_xyz > 0 .and. i_xyz <= this % n_cells) - end function valid_inds_rect + subroutine lattice_distance(this, xyz, uvw, i_xyz, d, lattice_trans) + class(Lattice), intent(in) :: this + real(C_DOUBLE), intent(in) :: xyz(3) + real(C_DOUBLE), intent(in) :: uvw(3) + integer(C_INT), intent(in) :: i_xyz(3) + real(C_DOUBLE), intent(out) :: d + integer(C_INT), intent(out) :: lattice_trans(3) + call lattice_distance_c(this % ptr, xyz, uvw, i_xyz, d, lattice_trans) + end subroutine lattice_distance + + function lattice_get(this, i_xyz) result(univ) + class(Lattice), intent(in) :: this + integer(C_INT), intent(in) :: i_xyz(3) + integer(C_INT32_T) :: univ + univ = lattice_universe_c(this % ptr, i_xyz) + end function lattice_get + + function lattice_get_indices(this, xyz) result(i_xyz) + class(Lattice), intent(in) :: this + real(C_DOUBLE), intent(in) :: xyz(3) + integer(C_INT) :: i_xyz(3) + call lattice_get_indices_c(this % ptr, xyz, i_xyz) + end function lattice_get_indices + + function lattice_get_local_xyz(this, global_xyz, i_xyz) & + result(local_xyz) + class(Lattice), intent(in) :: this + real(C_DOUBLE), intent(in) :: global_xyz(3) + integer(C_INT), intent(in) :: i_xyz(3) + real(C_DOUBLE) :: local_xyz(3) + call lattice_get_local_xyz_c(this % ptr, global_xyz, i_xyz, local_xyz) + end function lattice_get_local_xyz + + function lattice_offset(this, map, i_xyz) result(offset) + class(Lattice), intent(in) :: this + integer(C_INT), intent(in) :: map + integer(C_INT), intent(in) :: i_xyz(3) + integer(C_INT32_T) :: offset + offset = lattice_offset_c(this % ptr, map, i_xyz) + end function lattice_offset + + function lattice_outer(this) result(outer) + class(Lattice), intent(in) :: this + integer(C_INT32_T) :: outer + outer = lattice_outer_c(this % ptr) + end function lattice_outer + + subroutine lattice_to_hdf5(this, group) + class(Lattice), intent(in) :: this + integer(HID_T), intent(in) :: group + call lattice_to_hdf5_c(this % ptr, group) + end subroutine lattice_to_hdf5 !=============================================================================== - function valid_inds_hex(this, i_xyz) result(is_valid) - class(HexLattice), intent(in) :: this - integer, intent(in) :: i_xyz(3) - logical :: is_valid + function cell_id(this) result(id) + class(Cell), intent(in) :: this + integer(C_INT32_T) :: id + id = cell_id_c(this % ptr) + end function cell_id - is_valid = (all(i_xyz > 0) .and. & - &i_xyz(1) < 2*this % n_rings .and. & - &i_xyz(2) < 2*this % n_rings .and. & - &i_xyz(1) + i_xyz(2) > this % n_rings .and. & - &i_xyz(1) + i_xyz(2) < 3*this % n_rings .and. & - &i_xyz(3) <= this % n_axial) - end function valid_inds_hex + subroutine cell_set_id(this, id) + class(Cell), intent(in) :: this + integer(C_INT32_T), intent(in) :: id + call cell_set_id_c(this % ptr, id) + end subroutine cell_set_id -!=============================================================================== + function cell_type(this) result(type) + class(Cell), intent(in) :: this + integer(C_INT) :: type + type = cell_type_c(this % ptr) + end function cell_type - function get_inds_rect(this, global_xyz) result(i_xyz) - class(RectLattice), intent(in) :: this - real(8), intent(in) :: global_xyz(3) - integer :: i_xyz(3) + subroutine cell_set_type(this, type) + class(Cell), intent(in) :: this + integer(C_INT), intent(in) :: type + call cell_set_type_c(this % ptr, type) + end subroutine cell_set_type - real(8) :: xyz(3) ! global_xyz alias + function cell_universe(this) result(universe) + class(Cell), intent(in) :: this + integer(C_INT32_T) :: universe + universe = cell_universe_c(this % ptr) + end function cell_universe - xyz = global_xyz + subroutine cell_set_universe(this, universe) + class(Cell), intent(in) :: this + integer(C_INT32_T), intent(in) :: universe + call cell_set_universe_c(this % ptr, universe) + end subroutine cell_set_universe - i_xyz(1) = ceiling((xyz(1) - this % lower_left(1))/this % pitch(1)) - i_xyz(2) = ceiling((xyz(2) - this % lower_left(2))/this % pitch(2)) - if (this % is_3d) then - i_xyz(3) = ceiling((xyz(3) - this % lower_left(3))/this % pitch(3)) - else - i_xyz(3) = 1 - end if - end function get_inds_rect + function cell_fill(this) result(fill) + class(Cell), intent(in) :: this + integer(C_INT32_T) :: fill + fill = cell_fill_c(this % ptr) + end function cell_fill -!=============================================================================== + function cell_n_instances(this) result(n_instances) + class(Cell), intent(in) :: this + integer(C_INT32_T) :: n_instances + n_instances = cell_n_instances_c(this % ptr) + end function cell_n_instances - function get_inds_hex(this, global_xyz) result(i_xyz) - class(HexLattice), intent(in) :: this - real(8), intent(in) :: global_xyz(3) - integer :: i_xyz(3) + function cell_simple(this) result(simple) + class(Cell), intent(in) :: this + logical(C_BOOL) :: simple + simple = cell_simple_c(this % ptr) + end function cell_simple - real(8) :: xyz(3) ! global xyz relative to the center - real(8) :: alpha ! Skewed coord axis - real(8) :: xyz_t(3) ! Local xyz - real(8) :: d, d_min ! Squared distance from cell centers - integer :: i, j, k ! Iterators - integer :: k_min ! Minimum distance index + subroutine cell_distance(this, xyz, uvw, on_surface, min_dist, i_surf) + class(Cell), intent(in) :: this + real(C_DOUBLE), intent(in) :: xyz(3) + real(C_DOUBLE), intent(in) :: uvw(3) + integer(C_INT32_T), intent(in) :: on_surface + real(C_DOUBLE), intent(out) :: min_dist + integer(C_INT32_T), intent(out) :: i_surf + call cell_distance_c(this % ptr, xyz, uvw, on_surface, min_dist, i_surf) + end subroutine cell_distance - xyz(1) = global_xyz(1) - this % center(1) - xyz(2) = global_xyz(2) - this % center(2) + function cell_offset(this, map) result(offset) + class(Cell), intent(in) :: this + integer(C_INT), intent(in) :: map + integer(C_INT32_T) :: offset + offset = cell_offset_c(this % ptr, map) + end function cell_offset - ! Index z direction. - if (this % is_3d) then - xyz(3) = global_xyz(3) - this % center(3) - i_xyz(3) = ceiling(xyz(3)/this % pitch(2) + HALF*this % n_axial) - else - xyz(3) = global_xyz(3) - i_xyz(3) = 1 - end if - - ! Convert coordinates into skewed bases. The (x, alpha) basis is used to - ! find the index of the global coordinates to within 4 cells. - alpha = xyz(2) - xyz(1) / sqrt(THREE) - i_xyz(1) = floor(xyz(1) / (sqrt(THREE) / TWO * this % pitch(1))) - i_xyz(2) = floor(alpha / this % pitch(1)) - - ! Add offset to indices (the center cell is (i_x, i_alpha) = (0, 0) but - ! the array is offset so that the indices never go below 1). - i_xyz(1) = i_xyz(1) + this % n_rings - i_xyz(2) = i_xyz(2) + this % n_rings - - ! Calculate the (squared) distance between the particle and the centers of - ! the four possible cells. Regular hexagonal tiles form a centroidal - ! Voronoi tessellation so the global xyz should be in the hexagonal cell - ! that it is closest to the center of. This method is used over a - ! method that uses the remainders of the floor divisions above because it - ! provides better finite precision performance. Squared distances are - ! used becasue they are more computationally efficient than normal - ! distances. - k = 1 - d_min = INFINITY - do i = 0, 1 - do j = 0, 1 - xyz_t = this % get_local_xyz(global_xyz, i_xyz + [j, i, 0]) - d = xyz_t(1)**2 + xyz_t(2)**2 - if (d < d_min) then - d_min = d - k_min = k - end if - k = k + 1 - end do - end do - - ! Select the minimum squared distance which corresponds to the cell the - ! coordinates are in. - if (k_min == 2) then - i_xyz(1) = i_xyz(1) + 1 - else if (k_min == 3) then - i_xyz(2) = i_xyz(2) + 1 - else if (k_min == 4) then - i_xyz(1) = i_xyz(1) + 1 - i_xyz(2) = i_xyz(2) + 1 - end if - end function get_inds_hex - -!=============================================================================== - - function get_local_rect(this, global_xyz, i_xyz) result(local_xyz) - class(RectLattice), intent(in) :: this - real(8), intent(in) :: global_xyz(3) - integer, intent(in) :: i_xyz(3) - real(8) :: local_xyz(3) - - real(8) :: xyz(3) ! global_xyz alias - - xyz = global_xyz - - local_xyz(1) = xyz(1) - (this % lower_left(1) + & - (i_xyz(1) - HALF)*this % pitch(1)) - local_xyz(2) = xyz(2) - (this % lower_left(2) + & - (i_xyz(2) - HALF)*this % pitch(2)) - if (this % is_3d) then - local_xyz(3) = xyz(3) - (this % lower_left(3) + & - (i_xyz(3) - HALF)*this % pitch(3)) - else - local_xyz(3) = xyz(3) - end if - end function get_local_rect - -!=============================================================================== - - function get_local_hex(this, global_xyz, i_xyz) result(local_xyz) - class(HexLattice), intent(in) :: this - real(8), intent(in) :: global_xyz(3) - integer, intent(in) :: i_xyz(3) - real(8) :: local_xyz(3) - - real(8) :: xyz(3) ! global_xyz alias - - xyz = global_xyz - - ! x_l = x_g - (center + pitch_x*cos(30)*index_x) - local_xyz(1) = xyz(1) - (this % center(1) + & - sqrt(THREE) / TWO * (i_xyz(1) - this % n_rings) * this % pitch(1)) - ! y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y) - local_xyz(2) = xyz(2) - (this % center(2) + & - (i_xyz(2) - this % n_rings) * this % pitch(1) + & - (i_xyz(1) - this % n_rings) * this % pitch(1) / TWO) - if (this % is_3d) then - local_xyz(3) = xyz(3) - this % center(3) & - + (HALF*this % n_axial - i_xyz(3) + HALF) * this % pitch(2) - else - local_xyz(3) = xyz(3) - end if - end function get_local_hex + subroutine cell_to_hdf5(this, group) + class(Cell), intent(in) :: this + integer(HID_T), intent(in) :: group + call cell_to_hdf5_c(this % ptr, group) + end subroutine cell_to_hdf5 !=============================================================================== ! GET_TEMPERATURES returns a list of temperatures that each nuclide/S(a,b) table @@ -408,10 +516,15 @@ contains !=============================================================================== subroutine free_memory_geometry() + interface + subroutine free_memory_geometry_c() bind(C) + end subroutine free_memory_geometry_c + end interface + + call free_memory_geometry_c() n_cells = 0 n_universes = 0 - n_lattices = 0 if (allocated(cells)) deallocate(cells) if (allocated(universes)) deallocate(universes) @@ -432,6 +545,7 @@ contains 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_INT32_T) :: i integer(C_INT) :: err type(Cell), allocatable :: temp(:) ! temporary cells array @@ -453,7 +567,12 @@ contains ! Return indices in cells array if (present(index_start)) index_start = n_cells + 1 if (present(index_end)) index_end = n_cells + n - n_cells = n_cells + n + + ! Extend the C++ cells array and get pointers to the C++ objects + call extend_cells_c(n) + do i = n_cells - n, n_cells + cells(i) % ptr = cell_pointer_c(i - 1) + end do err = 0 end function openmc_extend_cells @@ -490,14 +609,14 @@ contains err = 0 if (index >= 1 .and. index <= size(cells)) then associate (c => cells(index)) - type = c % type + type = c % type() select case (type) case (FILL_MATERIAL) n = size(c % material) indices = C_LOC(c % material(1)) case (FILL_UNIVERSE, FILL_LATTICE) n = 1 - indices = C_LOC(c % fill) + indices = cell_fill_ptr(c % ptr) end select end associate else @@ -514,7 +633,7 @@ contains integer(C_INT) :: err if (index >= 1 .and. index <= size(cells)) then - id = cells(index) % id + id = cells(index) % id() err = 0 else err = E_OUT_OF_BOUNDS @@ -541,7 +660,7 @@ contains if (allocated(c % material)) deallocate(c % material) allocate(c % material(n)) - c % type = FILL_MATERIAL + call c % set_type(FILL_MATERIAL) do i = 1, n j = indices(i) if ((j >= 1 .and. j <= n_materials) .or. j == MATERIAL_VOID) then @@ -553,9 +672,9 @@ contains end if end do case (FILL_UNIVERSE) - c % type = FILL_UNIVERSE + call c % set_type(FILL_UNIVERSE) case (FILL_LATTICE) - c % type = FILL_LATTICE + call c % set_type(FILL_LATTICE) end select end associate else @@ -573,7 +692,7 @@ contains integer(C_INT) :: err if (index >= 1 .and. index <= n_cells) then - cells(index) % id = id + call cells(index) % set_id(id) call cell_dict % set(id, index) err = 0 else @@ -609,7 +728,7 @@ contains if (index >= 1 .and. index <= size(cells)) then ! error if the cell is filled with another universe - if (cells(index) % fill /= NONE) then + if (cells(index) % fill() /= C_NONE) then err = E_GEOMETRY call set_errmsg("Cannot set temperature on a cell filled & &with a universe.") diff --git a/src/hdf5_interface.h b/src/hdf5_interface.h index fd03eb4c67..59a914f52e 100644 --- a/src/hdf5_interface.h +++ b/src/hdf5_interface.h @@ -115,12 +115,21 @@ extern "C" void write_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_ const double* results); template void -write_double_1D(hid_t group_id, char const *name, - std::array &buffer) +write_int(hid_t group_id, char const *name, + const std::array &buffer, bool indep) +{ + hsize_t dims[1] {array_len}; + write_dataset(group_id, 1, dims, name, H5T_NATIVE_INT, buffer.data(), indep); +} + + +template void +write_double(hid_t group_id, char const *name, + const std::array &buffer, bool indep) { hsize_t dims[1] {array_len}; write_dataset(group_id, 1, dims, name, H5T_NATIVE_DOUBLE, - buffer.data(), false); + buffer.data(), indep); } } // namespace openmc diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 189ca7c212..6a9789b2e3 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -11,8 +11,7 @@ module input_xml use distribution_univariate use endf, only: reaction_name use error, only: fatal_error, warning, write_message, openmc_err_msg - use geometry, only: calc_offsets, maximum_levels, count_instance, & - neighbor_lists + use geometry, only: neighbor_lists use geometry_header use hdf5_interface use list_header, only: ListChar, ListInt, ListReal @@ -48,11 +47,50 @@ module input_xml save interface - subroutine read_surfaces(node_ptr) bind(C, name='read_surfaces') - use ISO_C_BINDING - implicit none + subroutine adjust_indices_c() bind(C) + end subroutine adjust_indices_c + + subroutine allocate_offset_tables(n_maps) bind(C) + import C_INT + integer(C_INT), intent(in), value :: n_maps + end subroutine allocate_offset_tables + + subroutine fill_offset_tables(target_univ_id, map) bind(C) + import C_INT32_T, C_INT + integer(C_INT32_T), intent(in), value :: target_univ_id + integer(C_INT), intent(in), value :: map + end subroutine fill_offset_tables + + subroutine count_cell_instances(univ_indx) bind(C) + import C_INT32_T + integer(C_INT32_T), intent(in), value :: univ_indx + end subroutine count_cell_instances + + subroutine read_surfaces(node_ptr) bind(C) + import C_PTR type(C_PTR) :: node_ptr end subroutine read_surfaces + + subroutine read_cells(node_ptr) bind(C) + import C_PTR + type(C_PTR) :: node_ptr + end subroutine read_cells + + subroutine read_lattices(node_ptr) bind(C) + import C_PTR + type(C_PTR) :: node_ptr + end subroutine read_lattices + + function find_root_universe() bind(C) result(root) + import C_INT32_T + integer(C_INT32_T) :: root + end function find_root_universe + + function maximum_levels(univ) bind(C) result(n) + import C_INT32_T, C_INT + integer(C_INT32_T), intent(in), value :: univ + integer(C_INT) :: n + end function maximum_levels end interface contains @@ -122,7 +160,7 @@ contains ! Perform some final operations to set up the geometry call adjust_indices() - call count_instance(universes(root_universe)) + call count_cell_instances(root_universe-1) ! After reading input and basic geometry setup is complete, build lists of ! neighboring cells for efficient tracking @@ -137,7 +175,7 @@ contains ! Check to make sure there are not too many nested coordinate levels in the ! geometry since the coordinate list is statically allocated for performance ! reasons - if (maximum_levels(universes(root_universe)) > MAX_COORD) then + if (maximum_levels(root_universe - 1) > MAX_COORD) then call fatal_error("Too many nested coordinate levels in the geometry. & &Try increasing the maximum number of coordinate levels by & &providing the CMake -Dmaxcoord= option.") @@ -913,12 +951,11 @@ contains subroutine read_geometry_xml() - integer :: i, j, k, m, i_x, i_a, input_index - integer :: n, n_mats, n_x, n_y, n_z, n_rings, n_rlats, n_hlats + integer :: i, j, k + integer :: n, n_mats, n_rlats, n_hlats integer :: id integer :: univ_id integer :: n_cells_in_univ - integer, allocatable :: temp_int_array(:) real(8) :: phi, theta, psi logical :: file_exists logical :: boundary_exists @@ -935,8 +972,6 @@ contains type(XMLNode), allocatable :: node_rlat_list(:) type(XMLNode), allocatable :: node_hlat_list(:) type(VectorInt) :: tokens - type(VectorInt) :: rpn - type(VectorInt) :: fill_univ_ids ! List of fill universe IDs type(VectorInt) :: univ_ids ! List of all universe IDs type(DictIntInt) :: cells_in_univ_dict ! Used to count how many cells each ! universe contains @@ -988,6 +1023,8 @@ contains ! ========================================================================== ! READ CELLS FROM GEOMETRY.XML + call read_cells(root % ptr) + ! Get pointer to list of XML call get_node_list(root, "cell", node_cell_list) @@ -1011,42 +1048,18 @@ contains do i = 1, n_cells c => cells(i) + c % ptr = cell_pointer_c(i - 1) + ! Initialize distribcell instances and distribcell index - c % instances = 0 c % distribcell_index = NONE ! Get pointer to i-th cell node node_cell = node_cell_list(i) - ! Copy data into cells - if (check_for_node(node_cell, "id")) then - call get_node_value(node_cell, "id", c % id) - else - call fatal_error("Must specify id of cell in geometry XML file.") - end if - - ! Copy cell name - if (check_for_node(node_cell, "name")) then - call get_node_value(node_cell, "name", c % name) - end if - - if (check_for_node(node_cell, "universe")) then - call get_node_value(node_cell, "universe", c % universe) - else - c % universe = 0 - end if - if (check_for_node(node_cell, "fill")) then - call get_node_value(node_cell, "fill", c % fill) - if (find(fill_univ_ids, c % fill) == -1) & - call fill_univ_ids % push_back(c % fill) - else - c % fill = NONE - end if - ! Check to make sure 'id' hasn't been used - if (cell_dict % has(c % id)) then + if (cell_dict % has(c % id())) then call fatal_error("Two or more cells use the same unique ID: " & - // to_str(c % id)) + // to_str(c % id())) end if ! Read material @@ -1068,7 +1081,7 @@ contains ! Check for error if (c % material(j) == ERROR_INT) then call fatal_error("Invalid material specified on cell " & - // to_str(c % id)) + // to_str(c % id())) end if end select end do @@ -1085,18 +1098,6 @@ contains c % material(1) = NONE end if - ! Check to make sure that either material or fill was specified - if (c % material(1) == NONE .and. c % fill == NONE) then - call fatal_error("Neither material nor fill was specified for cell " & - // trim(to_str(c % id))) - end if - - ! Check to make sure that both material and fill haven't been - ! specified simultaneously - if (c % material(1) /= NONE .and. c % fill /= NONE) then - call fatal_error("Cannot specify material and fill simultaneously") - end if - ! Check for region specification (also under deprecated name surfaces) if (check_for_node(node_cell, "surfaces")) then call warning("The use of 'surfaces' is deprecated and will be & @@ -1126,42 +1127,28 @@ contains end if end do - ! Use shunting-yard algorithm to determine RPN for surface algorithm - call generate_rpn(c%id, tokens, rpn) - ! Copy region spec and RPN form to cell arrays allocate(c % region(tokens%size())) - allocate(c % rpn(rpn%size())) c % region(:) = tokens%data(1:tokens%size()) - c % rpn(:) = rpn%data(1:rpn%size()) call tokens%clear() - call rpn%clear() end if if (.not. allocated(c%region)) allocate(c%region(0)) - if (.not. allocated(c%rpn)) allocate(c%rpn(0)) - - ! Check if this is a simple cell - if (any(c%rpn == OP_COMPLEMENT) .or. any(c%rpn == OP_UNION)) then - c%simple = .false. - else - c%simple = .true. - end if ! Rotation matrix if (check_for_node(node_cell, "rotation")) then ! Rotations can only be applied to cells that are being filled with ! another universe - if (c % fill == NONE) then + if (c % fill() == C_NONE) then call fatal_error("Cannot apply a rotation to cell " // trim(to_str(& - &c % id)) // " because it is not filled with another universe") + &c % id())) // " because it is not filled with another universe") end if ! Read number of rotation parameters n = node_word_count(node_cell, "rotation") if (n /= 3) then call fatal_error("Incorrect number of rotation parameters on cell " & - // to_str(c % id)) + // to_str(c % id())) end if ! Copy rotation angles in x,y,z directions @@ -1187,9 +1174,9 @@ contains if (check_for_node(node_cell, "translation")) then ! Translations can only be applied to cells that are being filled with ! another universe - if (c % fill == NONE) then + if (c % fill() == C_NONE) then call fatal_error("Cannot apply a translation to cell " & - // trim(to_str(c % id)) // " because it is not filled with & + // trim(to_str(c % id())) // " because it is not filled with & &another universe") end if @@ -1197,7 +1184,7 @@ contains n = node_word_count(node_cell, "translation") if (n /= 3) then call fatal_error("Incorrect number of translation parameters on & - &cell " // to_str(c % id)) + &cell " // to_str(c % id())) end if ! Copy translation vector @@ -1213,7 +1200,7 @@ contains if (n > 0) then ! Make sure this is a "normal" cell. if (c % material(1) == NONE) call fatal_error("Cell " & - // trim(to_str(c % id)) // " was specified with a temperature & + // trim(to_str(c % id())) // " was specified with a temperature & &but no material. Temperature specification is only valid for & &cells filled with a material.") @@ -1224,7 +1211,7 @@ contains ! Make sure all temperatues are positive do j = 1, size(c % sqrtkT) if (c % sqrtkT(j) < ZERO) call fatal_error("Cell " & - // trim(to_str(c % id)) // " was specified with a negative & + // trim(to_str(c % id())) // " was specified with a negative & &temperature. All cell temperatures must be non-negative.") end do @@ -1240,12 +1227,12 @@ contains end if ! Add cell to dictionary - call cell_dict % set(c % id, i) + call cell_dict % set(c % id(), i) ! For cells, we also need to check if there's a new universe -- ! also for every cell add 1 to the count of cells for the ! specified universe - univ_id = c % universe + univ_id = c % universe() if (.not. cells_in_univ_dict % has(univ_id)) then n_universes = n_universes + 1 n_cells_in_univ = 1 @@ -1261,6 +1248,8 @@ contains ! ========================================================================== ! READ LATTICES FROM GEOMETRY.XML + call read_lattices(root % ptr) + ! Get pointer to list of XML call get_node_list(root, "lattice", node_rlat_list) call get_node_list(root, "hex_lattice", node_hlat_list) @@ -1268,137 +1257,20 @@ contains ! Allocate lattices array n_rlats = size(node_rlat_list) n_hlats = size(node_hlat_list) - n_lattices = n_rlats + n_hlats - allocate(lattices(n_lattices)) + allocate(lattices(n_rlats + n_hlats)) RECT_LATTICES: do i = 1, n_rlats allocate(RectLattice::lattices(i) % obj) lat => lattices(i) % obj + lat % ptr = lattice_pointer_c(i - 1) select type(lat) type is (RectLattice) ! Get pointer to i-th lattice node_lat = node_rlat_list(i) - ! ID of lattice - if (check_for_node(node_lat, "id")) then - call get_node_value(node_lat, "id", lat % id) - else - call fatal_error("Must specify id of lattice in geometry XML file.") - end if - - ! Check to make sure 'id' hasn't been used - if (lattice_dict % has(lat % id)) then - call fatal_error("Two or more lattices use the same unique ID: " & - // to_str(lat % id)) - end if - - ! Copy lattice name - if (check_for_node(node_lat, "name")) then - call get_node_value(node_lat, "name", lat % name) - end if - - ! Read number of lattice cells in each dimension - n = node_word_count(node_lat, "dimension") - if (n == 2) then - call get_node_array(node_lat, "dimension", lat % n_cells(1:2)) - lat % n_cells(3) = 1 - lat % is_3d = .false. - else if (n == 3) then - call get_node_array(node_lat, "dimension", lat % n_cells) - lat % is_3d = .true. - else - call fatal_error("Rectangular lattice must be two or three dimensions.") - end if - - ! Read lattice lower-left location - if (node_word_count(node_lat, "lower_left") /= n) then - call fatal_error("Number of entries on must be the same & - &as the number of entries on .") - end if - - allocate(lat % lower_left(n)) - call get_node_array(node_lat, "lower_left", lat % lower_left) - - ! Read lattice pitches. - ! TODO: Remove this deprecation warning in a future release. - if (check_for_node(node_lat, "width")) then - call warning("The use of 'width' is deprecated and will be disallowed & - &in a future release. Use 'pitch' instead. The utility openmc/& - &src/utils/update_inputs.py can be used to automatically update & - &geometry.xml files.") - if (node_word_count(node_lat, "width") /= n) then - call fatal_error("Number of entries on must be the same as & - &the number of entries on .") - end if - - else if (node_word_count(node_lat, "pitch") /= n) then - call fatal_error("Number of entries on must be the same as & - &the number of entries on .") - end if - - allocate(lat % pitch(n)) - ! TODO: Remove the 'width' code in a future release. - if (check_for_node(node_lat, "width")) then - call get_node_array(node_lat, "width", lat % pitch) - else - call get_node_array(node_lat, "pitch", lat % pitch) - end if - - ! TODO: Remove deprecation warning in a future release. - if (check_for_node(node_lat, "type")) then - call warning("The use of 'type' is no longer needed. The utility & - &openmc/src/utils/update_inputs.py can be used to automatically & - &update geometry.xml files.") - end if - - ! Copy number of dimensions - n_x = lat % n_cells(1) - n_y = lat % n_cells(2) - n_z = lat % n_cells(3) - allocate(lat % universes(n_x, n_y, n_z)) - - ! Check that number of universes matches size - n = node_word_count(node_lat, "universes") - if (n /= n_x*n_y*n_z) then - call fatal_error("Number of universes on does not match & - &size of lattice " // trim(to_str(lat % id)) // ".") - end if - - allocate(temp_int_array(n)) - call get_node_array(node_lat, "universes", temp_int_array) - - ! Read universes - do m = 1, n_z - do k = 0, n_y - 1 - do j = 1, n_x - lat % universes(j, n_y - k, m) = & - temp_int_array(j + n_x*k + n_x*n_y*(m-1)) - if (find(fill_univ_ids, lat % universes(j, n_y - k, m)) == -1) & - call fill_univ_ids % push_back(lat % universes(j, n_y - k, m)) - end do - end do - end do - deallocate(temp_int_array) - - ! Read outer universe for area outside lattice. - lat % outer = NO_OUTER_UNIVERSE - if (check_for_node(node_lat, "outer")) then - call get_node_value(node_lat, "outer", lat % outer) - if (find(fill_univ_ids, lat % outer) == -1) & - call fill_univ_ids % push_back(lat % outer) - end if - - ! Check for 'outside' nodes which are no longer supported. - if (check_for_node(node_lat, "outside")) then - call fatal_error("The use of 'outside' in lattices is no longer & - &supported. Instead, use 'outer' which defines a universe rather & - &than a material. The utility openmc/src/utils/update_inputs.py & - &can be used automatically replace 'outside' with 'outer'.") - end if - ! Add lattice to dictionary - call lattice_dict % set(lat % id, i) + call lattice_dict % set(lat % id(), i) end select end do RECT_LATTICES @@ -1406,186 +1278,15 @@ contains HEX_LATTICES: do i = 1, n_hlats allocate(HexLattice::lattices(n_rlats + i) % obj) lat => lattices(n_rlats + i) % obj + lat % ptr = lattice_pointer_c(n_rlats + i - 1) select type (lat) type is (HexLattice) ! Get pointer to i-th lattice node_lat = node_hlat_list(i) - ! ID of lattice - if (check_for_node(node_lat, "id")) then - call get_node_value(node_lat, "id", lat % id) - else - call fatal_error("Must specify id of lattice in geometry XML file.") - end if - - ! Check to make sure 'id' hasn't been used - if (lattice_dict % has(lat % id)) then - call fatal_error("Two or more lattices use the same unique ID: " & - // to_str(lat % id)) - end if - - ! Copy lattice name - if (check_for_node(node_lat, "name")) then - call get_node_value(node_lat, "name", lat % name) - end if - - ! Read number of lattice cells in each dimension - call get_node_value(node_lat, "n_rings", lat % n_rings) - if (check_for_node(node_lat, "n_axial")) then - call get_node_value(node_lat, "n_axial", lat % n_axial) - lat % is_3d = .true. - else - lat % n_axial = 1 - lat % is_3d = .false. - end if - - ! Read lattice lower-left location - n = node_word_count(node_lat, "center") - if (lat % is_3d .and. n /= 3) then - call fatal_error("A hexagonal lattice with must have & - &
specified by 3 numbers.") - else if ((.not. lat % is_3d) .and. n /= 2) then - call fatal_error("A hexagonal lattice without must have & - &
specified by 2 numbers.") - end if - - allocate(lat % center(n)) - call get_node_array(node_lat, "center", lat % center) - - ! Read lattice pitches - n = node_word_count(node_lat, "pitch") - if (lat % is_3d .and. n /= 2) then - call fatal_error("A hexagonal lattice with must have & - &specified by 2 numbers.") - else if ((.not. lat % is_3d) .and. n /= 1) then - call fatal_error("A hexagonal lattice without must have & - & specified by 1 number.") - end if - - allocate(lat % pitch(n)) - call get_node_array(node_lat, "pitch", lat % pitch) - - ! Copy number of dimensions - n_rings = lat % n_rings - n_z = lat % n_axial - allocate(lat % universes(2*n_rings - 1, 2*n_rings - 1, n_z)) - - ! Check that number of universes matches size - n = node_word_count(node_lat, "universes") - if (n /= (3*n_rings**2 - 3*n_rings + 1)*n_z) then - call fatal_error("Number of universes on does not match & - &size of lattice " // trim(to_str(lat % id)) // ".") - end if - - allocate(temp_int_array(n)) - call get_node_array(node_lat, "universes", temp_int_array) - - ! Read universes - ! Universes in hexagonal lattices are stored in a manner that represents - ! a skewed coordinate system: (x, alpha) rather than (x, y). There is - ! no obvious, direct relationship between the order of universes in the - ! input and the order that they will be stored in the skewed array so - ! the following code walks a set of index values across the skewed array - ! in a manner that matches the input order. Note that i_x = 0, i_a = 0 - ! corresponds to the center of the hexagonal lattice. - - input_index = 1 - do m = 1, n_z - ! Initialize lattice indecies. - i_x = 1 - i_a = n_rings - 1 - - ! Map upper triangular region of hexagonal lattice. - do k = 1, n_rings-1 - ! Walk index to lower-left neighbor of last row start. - i_x = i_x - 1 - do j = 1, k - ! Place universe in array. - lat % universes(i_x + n_rings, i_a + n_rings, m) = & - temp_int_array(input_index) - if (find(fill_univ_ids, temp_int_array(input_index)) == -1) & - call fill_univ_ids % push_back(temp_int_array(input_index)) - ! Walk index to closest non-adjacent right neighbor. - i_x = i_x + 2 - i_a = i_a - 1 - ! Increment XML array index. - input_index = input_index + 1 - end do - ! Return lattice index to start of current row. - i_x = i_x - 2*k - i_a = i_a + k - end do - - ! Map middle square region of hexagonal lattice. - do k = 1, 2*n_rings - 1 - if (mod(k, 2) == 1) then - ! Walk index to lower-left neighbor of last row start. - i_x = i_x - 1 - else - ! Walk index to lower-right neighbor of last row start - i_x = i_x + 1 - i_a = i_a - 1 - end if - do j = 1, n_rings - mod(k-1, 2) - ! Place universe in array. - lat % universes(i_x + n_rings, i_a + n_rings, m) = & - temp_int_array(input_index) - if (find(fill_univ_ids, temp_int_array(input_index)) == -1) & - call fill_univ_ids % push_back(temp_int_array(input_index)) - ! Walk index to closest non-adjacent right neighbor. - i_x = i_x + 2 - i_a = i_a - 1 - ! Increment XML array index. - input_index = input_index + 1 - end do - ! Return lattice index to start of current row. - i_x = i_x - 2*(n_rings - mod(k-1, 2)) - i_a = i_a + n_rings - mod(k-1, 2) - end do - - ! Map lower triangular region of hexagonal lattice. - do k = 1, n_rings-1 - ! Walk index to lower-right neighbor of last row start. - i_x = i_x + 1 - i_a = i_a - 1 - do j = 1, n_rings - k - ! Place universe in array. - lat % universes(i_x + n_rings, i_a + n_rings, m) = & - temp_int_array(input_index) - if (find(fill_univ_ids, temp_int_array(input_index)) == -1) & - call fill_univ_ids % push_back(temp_int_array(input_index)) - ! Walk index to closest non-adjacent right neighbor. - i_x = i_x + 2 - i_a = i_a - 1 - ! Increment XML array index. - input_index = input_index + 1 - end do - ! Return lattice index to start of current row. - i_x = i_x - 2*(n_rings - k) - i_a = i_a + n_rings - k - end do - end do - deallocate(temp_int_array) - - ! Read outer universe for area outside lattice. - lat % outer = NO_OUTER_UNIVERSE - if (check_for_node(node_lat, "outer")) then - call get_node_value(node_lat, "outer", lat % outer) - if (find(fill_univ_ids, lat % outer) == -1) & - call fill_univ_ids % push_back(lat % outer) - end if - - ! Check for 'outside' nodes which are no longer supported. - if (check_for_node(node_lat, "outside")) then - call fatal_error("The use of 'outside' in lattices is no longer & - &supported. Instead, use 'outer' which defines a universe rather & - &than a material. The utility openmc/src/utils/update_inputs.py & - &can be used automatically replace 'outside' with 'outer'.") - end if - ! Add lattice to dictionary - call lattice_dict % set(lat % id, n_rlats + i) + call lattice_dict % set(lat % id(), n_rlats + i) end select end do HEX_LATTICES @@ -1603,23 +1304,13 @@ contains n_cells_in_univ = cells_in_univ_dict % get(u % id) allocate(u % cells(n_cells_in_univ)) u % cells(:) = 0 - - ! Check whether universe is a fill universe - if (find(fill_univ_ids, u % id) == -1) then - if (root_universe > 0) then - call fatal_error("Two or more universes are not used as fill & - &universes, so it is not possible to distinguish which one & - &is the root universe.") - else - root_universe = i - end if - end if end associate end do + root_universe = find_root_universe() + 1 do i = 1, n_cells ! Get index in universes array - j = universe_dict % get(cells(i) % universe) + j = universe_dict % get(cells(i) % universe()) ! Set the first zero entry in the universe cells array to the index in the ! global cells array @@ -3743,93 +3434,6 @@ contains end subroutine read_mg_cross_sections_header -!=============================================================================== -! GENERATE_RPN implements the shunting-yard algorithm to generate a Reverse -! Polish notation (RPN) expression for the region specification of a cell given -! the infix notation. -!=============================================================================== - - subroutine generate_rpn(cell_id, tokens, output) - integer, intent(in) :: cell_id - type(VectorInt), intent(in) :: tokens ! infix notation - type(VectorInt), intent(inout) :: output ! RPN notation - - integer :: i - integer :: token - integer :: op - type(VectorInt) :: stack - - do i = 1, tokens%size() - token = tokens%data(i) - - if (token < OP_UNION) then - ! If token is not an operator, add it to output - call output%push_back(token) - - elseif (token < OP_RIGHT_PAREN) then - ! Regular operators union, intersection, complement - do while (stack%size() > 0) - op = stack%data(stack%size()) - - if (op < OP_RIGHT_PAREN .and. & - ((token == OP_COMPLEMENT .and. token < op) .or. & - (token /= OP_COMPLEMENT .and. token <= op))) then - ! While there is an operator, op, on top of the stack, if the token - ! is left-associative and its precedence is less than or equal to - ! that of op or if the token is right-associative and its precedence - ! is less than that of op, move op to the output queue and push the - ! token on to the stack. Note that only complement is - ! right-associative. - call output%push_back(op) - call stack%pop_back() - else - exit - end if - end do - - call stack%push_back(token) - - elseif (token == OP_LEFT_PAREN) then - ! If the token is a left parenthesis, push it onto the stack - call stack%push_back(token) - - else - ! If the token is a right parenthesis, move operators from the stack to - ! the output queue until reaching the left parenthesis. - do - ! If we run out of operators without finding a left parenthesis, it - ! means there are mismatched parentheses. - if (stack%size() == 0) then - call fatal_error('Mimatched parentheses in region specification & - &for cell ' // trim(to_str(cell_id)) // '.') - end if - - op = stack%data(stack%size()) - if (op == OP_LEFT_PAREN) exit - call output%push_back(op) - call stack%pop_back() - end do - - ! Pop the left parenthesis. - call stack%pop_back() - end if - end do - - ! While there are operators on the stack, move them to the output queue - do while (stack%size() > 0) - op = stack%data(stack%size()) - - ! If the operator is a parenthesis, it is mismatched - if (op >= OP_RIGHT_PAREN) then - call fatal_error('Mimatched parentheses in region specification & - &for cell ' // trim(to_str(cell_id)) // '.') - end if - - call output%push_back(op) - call stack%pop_back() - end do - end subroutine generate_rpn - !=============================================================================== ! NORMALIZE_AO Normalize the nuclide atom percents !=============================================================================== @@ -4150,117 +3754,31 @@ contains integer :: i ! index for various purposes integer :: j ! index for various purposes - integer :: k ! loop index for lattices - integer :: m ! loop index for lattices - integer :: lid ! lattice IDs integer :: id ! user-specified id - class(Lattice), pointer :: lat => null() + + call adjust_indices_c() do i = 1, n_cells - ! ======================================================================= - ! ADJUST UNIVERSE INDEX FOR EACH CELL associate (c => cells(i)) - id = c % universe - if (universe_dict % has(id)) then - c % universe = universe_dict % get(id) - else - call fatal_error("Could not find universe " // trim(to_str(id)) & - &// " specified on cell " // trim(to_str(c % id))) - end if - ! ======================================================================= ! ADJUST MATERIAL/FILL POINTERS FOR EACH CELL - if (c % material(1) == NONE) then - id = c % fill - if (universe_dict % has(id)) then - c % type = FILL_UNIVERSE - c % fill = universe_dict % get(id) - elseif (lattice_dict % has(id)) then - lid = lattice_dict % get(id) - c % type = FILL_LATTICE - c % fill = lid - else - call fatal_error("Specified fill " // trim(to_str(id)) // " on cell "& - // trim(to_str(c % id)) // " is neither a universe nor a & - &lattice.") - end if - else + if (c % material(1) /= NONE) then do j = 1, size(c % material) id = c % material(j) if (id == MATERIAL_VOID) then - c % type = FILL_MATERIAL else if (material_dict % has(id)) then - c % type = FILL_MATERIAL c % material(j) = material_dict % get(id) else call fatal_error("Could not find material " // trim(to_str(id)) & - // " specified on cell " // trim(to_str(c % id))) + // " specified on cell " // trim(to_str(c % id()))) end if end do end if end associate end do - ! ========================================================================== - ! ADJUST UNIVERSE INDICES FOR EACH LATTICE - - do i = 1, n_lattices - lat => lattices(i) % obj - select type (lat) - - type is (RectLattice) - do m = 1, lat % n_cells(3) - do k = 1, lat % n_cells(2) - do j = 1, lat % n_cells(1) - id = lat % universes(j,k,m) - if (universe_dict % has(id)) then - lat % universes(j,k,m) = universe_dict % get(id) - else - call fatal_error("Invalid universe number " & - &// trim(to_str(id)) // " specified on lattice " & - &// trim(to_str(lat % id))) - end if - end do - end do - end do - - type is (HexLattice) - do m = 1, lat % n_axial - do k = 1, 2*lat % n_rings - 1 - do j = 1, 2*lat % n_rings - 1 - if (j + k < lat % n_rings + 1) then - cycle - else if (j + k > 3*lat % n_rings - 1) then - cycle - end if - id = lat % universes(j, k, m) - if (universe_dict % has(id)) then - lat % universes(j, k, m) = universe_dict % get(id) - else - call fatal_error("Invalid universe number " & - &// trim(to_str(id)) // " specified on lattice " & - &// trim(to_str(lat % id))) - end if - end do - end do - end do - - end select - - if (lat % outer /= NO_OUTER_UNIVERSE) then - if (universe_dict % has(lat % outer)) then - lat % outer = universe_dict % get(lat % outer) - else - call fatal_error("Invalid universe number " & - &// trim(to_str(lat % outer)) & - &// " specified on lattice " // trim(to_str(lat % id))) - end if - end if - - end do - end subroutine adjust_indices !=============================================================================== @@ -4270,94 +3788,11 @@ contains subroutine prepare_distribcell() - integer :: i, j ! Tally, filter loop counters - logical :: distribcell_active ! Does simulation use distribcell? - integer, allocatable :: univ_list(:) ! Target offsets - integer, allocatable :: counts(:,:) ! Target count - logical, allocatable :: found(:,:) ! Target found + integer :: i, j, k + type(SetInt) :: cell_list ! distribcells to track + type(ListInt) :: univ_list ! universes containing distribcells - ! Assume distribcell is not needed until proven otherwise. - distribcell_active = .false. - - ! We need distribcell if any tallies have distribcell filters. - do i = 1, n_tallies - do j = 1, size(tallies(i) % obj % filter) - select type(filt => filters(tallies(i) % obj % filter(j)) % obj) - type is (DistribcellFilter) - distribcell_active = .true. - end select - end do - end do - - ! We also need distribcell if any distributed materials or distributed - ! temperatues are present. - if (.not. distribcell_active) then - do i = 1, n_cells - if (size(cells(i) % material) > 1 .or. size(cells(i) % sqrtkT) > 1) then - distribcell_active = .true. - exit - end if - end do - end if - - ! If distribcell isn't used in this simulation then no more work left to do. - if (.not. distribcell_active) return - - ! Make sure the number of materials and temperatures matches the number of - ! cell instances. - do i = 1, n_cells - associate (c => cells(i)) - if (size(c % material) > 1) then - if (size(c % material) /= c % instances) then - call fatal_error("Cell " // trim(to_str(c % id)) // " was & - &specified with " // trim(to_str(size(c % material))) & - // " materials but has " // trim(to_str(c % instances)) & - // " distributed instances. The number of materials must & - &equal one or the number of instances.") - end if - end if - if (size(c % sqrtkT) > 1) then - if (size(c % sqrtkT) /= c % instances) then - call fatal_error("Cell " // trim(to_str(c % id)) // " was & - &specified with " // trim(to_str(size(c % sqrtkT))) & - // " temperatures but has " // trim(to_str(c % instances)) & - // " distributed instances. The number of temperatures must & - &equal one or the number of instances.") - end if - end if - end associate - end do - - ! Allocate offset maps at each level in the geometry - call allocate_offsets(univ_list, counts, found) - - ! Calculate offsets for each target distribcell - do i = 1, n_maps - do j = 1, n_universes - call calc_offsets(univ_list(i), i, universes(j), counts, found) - end do - end do - - end subroutine prepare_distribcell - -!=============================================================================== -! ALLOCATE_OFFSETS determines the number of maps needed and allocates required -! memory for distribcell offset tables -!=============================================================================== - - recursive subroutine allocate_offsets(univ_list, counts, found) - - integer, intent(out), allocatable :: univ_list(:) ! Target offsets - integer, intent(out), allocatable :: counts(:,:) ! Target count - logical, intent(out), allocatable :: found(:,:) ! Target found - - integer :: i, j, k ! Loop counters - type(SetInt) :: cell_list ! distribells to track - - ! Begin gathering list of cells in distribcell tallies - n_maps = 0 - - ! List all cells referenced in distribcell filters. + ! Find all cells listed in a distribcell filter. do i = 1, n_tallies do j = 1, size(tallies(i) % obj % filter) select type(filt => filters(tallies(i) % obj % filter(j)) % obj) @@ -4367,35 +3802,38 @@ contains end do end do - ! List all cells with multiple (distributed) materials or temperatures. + ! Find all cells with multiple (distributed) materials or temperatures. do i = 1, n_cells if (size(cells(i) % material) > 1 .or. size(cells(i) % sqrtkT) > 1) then call cell_list % add(i) end if end do - ! Compute the number of unique universes containing these distribcells - ! to determine the number of offset tables to allocate - do i = 1, n_universes - do j = 1, size(universes(i) % cells) - if (cell_list % contains(universes(i) % cells(j))) then - n_maps = n_maps + 1 + ! Make sure the number of distributed materials and temperatures matches the + ! number of respective cell instances. + do i = 1, n_cells + associate (c => cells(i)) + if (size(c % material) > 1) then + if (size(c % material) /= c % n_instances()) then + call fatal_error("Cell " // trim(to_str(c % id())) // " was & + &specified with " // trim(to_str(size(c % material))) & + // " materials but has " // trim(to_str(c % n_instances())) & + // " distributed instances. The number of materials must & + &equal one or the number of instances.") + end if end if - end do + if (size(c % sqrtkT) > 1) then + if (size(c % sqrtkT) /= c % n_instances()) then + call fatal_error("Cell " // trim(to_str(c % id())) // " was & + &specified with " // trim(to_str(size(c % sqrtkT))) & + // " temperatures but has " // trim(to_str(c % n_instances()))& + // " distributed instances. The number of temperatures must & + &equal one or the number of instances.") + end if + end if + end associate end do - ! Allocate the list of offset tables for each unique universe - allocate(univ_list(n_maps)) - - ! Allocate list to accumulate target distribcell counts in each universe - allocate(counts(n_universes, n_maps)) - counts(:,:) = 0 - - ! Allocate list to track if target distribcells are found in each universe - allocate(found(n_universes, n_maps)) - found(:,:) = .false. - - ! Search through universes for distributed cells and assign each one a ! unique distribcell array index. k = 1 @@ -4403,39 +3841,18 @@ contains do j = 1, size(universes(i) % cells) if (cell_list % contains(universes(i) % cells(j))) then cells(universes(i) % cells(j)) % distribcell_index = k - univ_list(k) = universes(i) % id + call univ_list % append(universes(i) % id) k = k + 1 end if end do end do - ! Allocate the offset tables for lattices - do i = 1, n_lattices - associate(lat => lattices(i) % obj) - select type(lat) - - type is (RectLattice) - allocate(lat % offset(n_maps, lat % n_cells(1), lat % n_cells(2), & - lat % n_cells(3))) - type is (HexLattice) - allocate(lat % offset(n_maps, 2 * lat % n_rings - 1, & - 2 * lat % n_rings - 1, lat % n_axial)) - end select - - lat % offset(:, :, :, :) = 0 - end associate + ! Allocate and fill cell and lattice offset tables. + call allocate_offset_tables(univ_list % size()) + do i = 1, univ_list % size() + call fill_offset_tables(univ_list % get_item(i), i-1) end do - ! Allocate offset table for fill cells - do i = 1, n_cells - if (cells(i) % type /= FILL_MATERIAL) then - allocate(cells(i) % offset(n_maps)) - end if - end do - - ! Free up memory - call cell_list % clear() - - end subroutine allocate_offsets + end subroutine prepare_distribcell end module input_xml diff --git a/src/lattice.cpp b/src/lattice.cpp new file mode 100644 index 0000000000..959ab932e4 --- /dev/null +++ b/src/lattice.cpp @@ -0,0 +1,946 @@ +#include "lattice.h" + +#include +#include +#include + +#include "cell.h" +#include "error.h" +#include "geometry_aux.h" +#include "hdf5_interface.h" +#include "string_utils.h" +#include "xml_interface.h" + + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +std::vector lattices_c; + +std::unordered_map lattice_map; + +//============================================================================== +// Lattice implementation +//============================================================================== + +Lattice::Lattice(pugi::xml_node lat_node) +{ + if (check_for_node(lat_node, "id")) { + id = stoi(get_node_value(lat_node, "id")); + } else { + fatal_error("Must specify id of lattice in geometry XML file."); + } + + if (check_for_node(lat_node, "name")) { + name = get_node_value(lat_node, "name"); + } + + if (check_for_node(lat_node, "outer")) { + outer = stoi(get_node_value(lat_node, "outer")); + } +} + +//============================================================================== + +LatticeIter Lattice::begin() +{return LatticeIter(*this, 0);} + +LatticeIter Lattice::end() +{return LatticeIter(*this, universes.size());} + +ReverseLatticeIter Lattice::rbegin() +{return ReverseLatticeIter(*this, universes.size()-1);} + +ReverseLatticeIter Lattice::rend() +{return ReverseLatticeIter(*this, -1);} + +//============================================================================== + +void +Lattice::adjust_indices() +{ + // Adjust the indices for the universes array. + for (LatticeIter it = begin(); it != end(); ++it) { + int uid = *it; + auto search = universe_map.find(uid); + if (search != universe_map.end()) { + *it = search->second; + } else { + std::stringstream err_msg; + err_msg << "Invalid universe number " << uid << " specified on " + "lattice " << id; + fatal_error(err_msg); + } + } + + // Adjust the index for the outer universe. + if (outer != NO_OUTER_UNIVERSE) { + auto search = universe_map.find(outer); + if (search != universe_map.end()) { + outer = search->second; + } else { + std::stringstream err_msg; + err_msg << "Invalid universe number " << outer << " specified on " + "lattice " << id; + fatal_error(err_msg); + } + } +} + +//============================================================================== + +int32_t +Lattice::fill_offset_table(int32_t offset, int32_t target_univ_id, int map) +{ + for (LatticeIter it = begin(); it != end(); ++it) { + offsets[map * universes.size() + it.indx] = offset; + offset += count_universe_instances(*it, target_univ_id); + } + return offset; +} + +//============================================================================== + +void +Lattice::to_hdf5(hid_t lattices_group) const +{ + // Make a group for the lattice. + std::string group_name {"lattice "}; + group_name += std::to_string(id); + hid_t lat_group = create_group(lattices_group, group_name); + + // Write the name and outer universe. + if (!name.empty()) { + write_string(lat_group, "name", name, false); + } + + if (outer != NO_OUTER_UNIVERSE) { + int32_t outer_id = global_universes[outer]->id; + write_int(lat_group, 0, nullptr, "outer", &outer_id, false); + } else { + write_int(lat_group, 0, nullptr, "outer", &outer, false); + } + + // Call subclass-overriden function to fill in other details. + to_hdf5_inner(lat_group); + + close_group(lat_group); +} + +//============================================================================== +// RectLattice implementation +//============================================================================== + +RectLattice::RectLattice(pugi::xml_node lat_node) + : Lattice {lat_node} +{ + // Read the number of lattice cells in each dimension. + std::string dimension_str {get_node_value(lat_node, "dimension")}; + std::vector dimension_words {split(dimension_str)}; + if (dimension_words.size() == 2) { + n_cells[0] = stoi(dimension_words[0]); + n_cells[1] = stoi(dimension_words[1]); + n_cells[2] = 1; + is_3d = false; + } else if (dimension_words.size() == 3) { + n_cells[0] = stoi(dimension_words[0]); + n_cells[1] = stoi(dimension_words[1]); + n_cells[2] = stoi(dimension_words[2]); + is_3d = true; + } else { + fatal_error("Rectangular lattice must be two or three dimensions."); + } + + // Read the lattice lower-left location. + std::string ll_str {get_node_value(lat_node, "lower_left")}; + std::vector ll_words {split(ll_str)}; + if (ll_words.size() != dimension_words.size()) { + fatal_error("Number of entries on must be the same as the " + "number of entries on ."); + } + lower_left[0] = stod(ll_words[0]); + lower_left[1] = stod(ll_words[1]); + if (is_3d) {lower_left[2] = stod(ll_words[2]);} + + // Read the lattice pitches. + std::string pitch_str {get_node_value(lat_node, "pitch")}; + std::vector pitch_words {split(pitch_str)}; + if (pitch_words.size() != dimension_words.size()) { + fatal_error("Number of entries on must be the same as the " + "number of entries on ."); + } + pitch[0] = stod(pitch_words[0]); + pitch[1] = stod(pitch_words[1]); + if (is_3d) {pitch[2] = stod(pitch_words[2]);} + + // Read the universes and make sure the correct number was specified. + std::string univ_str {get_node_value(lat_node, "universes")}; + std::vector univ_words {split(univ_str)}; + if (univ_words.size() != nx*ny*nz) { + std::stringstream err_msg; + err_msg << "Expected " << nx*ny*nz + << " universes for a rectangular lattice of size " + << nx << "x" << ny << "x" << nz << " but " << univ_words.size() + << " were specified."; + fatal_error(err_msg); + } + + // Parse the universes. + universes.resize(nx*ny*nz, C_NONE); + for (int iz = 0; iz < nz; iz++) { + for (int iy = ny-1; iy > -1; iy--) { + for (int ix = 0; ix < nx; ix++) { + int indx1 = nx*ny*iz + nx*(ny-iy-1) + ix; + int indx2 = nx*ny*iz + nx*iy + ix; + universes[indx1] = stoi(univ_words[indx2]); + } + } + } +} + +//============================================================================== + +int32_t& +RectLattice::operator[](const int i_xyz[3]) +{ + int indx = nx*ny*i_xyz[2] + nx*i_xyz[1] + i_xyz[0]; + return universes[indx]; +} + +//============================================================================== + +bool +RectLattice::are_valid_indices(const int i_xyz[3]) const +{ + return ( (i_xyz[0] >= 0) && (i_xyz[0] < n_cells[0]) + && (i_xyz[1] >= 0) && (i_xyz[1] < n_cells[1]) + && (i_xyz[2] >= 0) && (i_xyz[2] < n_cells[2])); +} + +//============================================================================== + +std::pair> +RectLattice::distance(const double xyz[3], const double uvw[3], + const int i_xyz[3]) const +{ + // Get short aliases to the coordinates. + double x {xyz[0]}; + double y {xyz[1]}; + double z {xyz[2]}; + double u {uvw[0]}; + double v {uvw[1]}; + + // Determine the oncoming edge. + double x0 {copysign(0.5 * pitch[0], u)}; + double y0 {copysign(0.5 * pitch[1], v)}; + + // Left and right sides + double d {INFTY}; + std::array lattice_trans; + if ((std::abs(x - x0) > FP_PRECISION) && u != 0) { + d = (x0 - x) / u; + if (u > 0) { + lattice_trans = {1, 0, 0}; + } else { + lattice_trans = {-1, 0, 0}; + } + } + + // Front and back sides + if ((std::abs(y - y0) > FP_PRECISION) && v != 0) { + double this_d = (y0 - y) / v; + if (this_d < d) { + d = this_d; + if (v > 0) { + lattice_trans = {0, 1, 0}; + } else { + lattice_trans = {0, -1, 0}; + } + } + } + + // Top and bottom sides + if (is_3d) { + double w {uvw[2]}; + double z0 {copysign(0.5 * pitch[2], w)}; + if ((std::abs(z - z0) > FP_PRECISION) && w != 0) { + double this_d = (z0 - z) / w; + if (this_d < d) { + d = this_d; + if (w > 0) { + lattice_trans = {0, 0, 1}; + } else { + lattice_trans = {0, 0, -1}; + } + } + } + } + + return {d, lattice_trans}; +} + +//============================================================================== + +std::array +RectLattice::get_indices(const double xyz[3]) const +{ + int ix {static_cast(std::ceil((xyz[0] - lower_left[0]) / pitch[0]))-1}; + int iy {static_cast(std::ceil((xyz[1] - lower_left[1]) / pitch[1]))-1}; + int iz; + if (is_3d) { + iz = static_cast(std::ceil((xyz[2] - lower_left[2]) / pitch[2]))-1; + } else { + iz = 0; + } + return {ix, iy, iz}; +} + +//============================================================================== + +std::array +RectLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const +{ + std::array local_xyz; + local_xyz[0] = global_xyz[0] - (lower_left[0] + (i_xyz[0] + 0.5)*pitch[0]); + local_xyz[1] = global_xyz[1] - (lower_left[1] + (i_xyz[1] + 0.5)*pitch[1]); + if (is_3d) { + local_xyz[2] = global_xyz[2] - (lower_left[2] + (i_xyz[2] + 0.5)*pitch[2]); + } else { + local_xyz[2] = global_xyz[2]; + } + return local_xyz; +} + +//============================================================================== + +int32_t& +RectLattice::offset(int map, const int i_xyz[3]) +{ + return offsets[nx*ny*nz*map + nx*ny*i_xyz[2] + nx*i_xyz[1] + i_xyz[0]]; +} + +//============================================================================== + +std::string +RectLattice::index_to_string(int indx) const +{ + int iz {indx / (nx * ny)}; + int iy {(indx - nx * ny * iz) / nx}; + int ix {indx - nx * ny * iz - nx * iy}; + std::string out {std::to_string(ix)}; + out += ','; + out += std::to_string(iy); + if (is_3d) { + out += ','; + out += std::to_string(iz); + } + return out; +} + +//============================================================================== + +void +RectLattice::to_hdf5_inner(hid_t lat_group) const +{ + // Write basic lattice information. + write_string(lat_group, "type", "rectangular", false); + if (is_3d) { + write_double(lat_group, "pitch", pitch, false); + write_double(lat_group, "lower_left", lower_left, false); + write_int(lat_group, "dimension", n_cells, false); + } else { + std::array pitch_short {{pitch[0], pitch[1]}}; + write_double(lat_group, "pitch", pitch_short, false); + std::array ll_short {{lower_left[0], lower_left[1]}}; + write_double(lat_group, "lower_left", ll_short, false); + std::array nc_short {{n_cells[0], n_cells[1]}}; + write_int(lat_group, "dimension", nc_short, false); + } + + // Write the universe ids. The convention here is to switch the ordering on + // the y-axis to match the way universes are input in a text file. + if (is_3d) { + hsize_t nx {static_cast(n_cells[0])}; + hsize_t ny {static_cast(n_cells[1])}; + hsize_t nz {static_cast(n_cells[2])}; + int out[nx*ny*nz]; + + for (int m = 0; m < nz; m++) { + for (int k = 0; k < ny; k++) { + for (int j = 0; j < nx; j++) { + int indx1 = nx*ny*m + nx*k + j; + int indx2 = nx*ny*m + nx*(ny-k-1) + j; + out[indx2] = global_universes[universes[indx1]]->id; + } + } + } + + hsize_t dims[3] {nz, ny, nx}; + write_int(lat_group, 3, dims, "universes", out, false); + + } else { + hsize_t nx {static_cast(n_cells[0])}; + hsize_t ny {static_cast(n_cells[1])}; + int out[nx*ny]; + + for (int k = 0; k < ny; k++) { + for (int j = 0; j < nx; j++) { + int indx1 = nx*k + j; + int indx2 = nx*(ny-k-1) + j; + out[indx2] = global_universes[universes[indx1]]->id; + } + } + + hsize_t dims[3] {1, ny, nx}; + write_int(lat_group, 3, dims, "universes", out, false); + } +} + +//============================================================================== +// HexLattice implementation +//============================================================================== + +HexLattice::HexLattice(pugi::xml_node lat_node) + : Lattice {lat_node} +{ + // Read the number of lattice cells in each dimension. + n_rings = stoi(get_node_value(lat_node, "n_rings")); + if (check_for_node(lat_node, "n_axial")) { + n_axial = stoi(get_node_value(lat_node, "n_axial")); + is_3d = true; + } else { + n_axial = 1; + is_3d = false; + } + + // Read the lattice center. + std::string center_str {get_node_value(lat_node, "center")}; + std::vector center_words {split(center_str)}; + if (is_3d && (center_words.size() != 3)) { + fatal_error("A hexagonal lattice with must have
" + "specified by 3 numbers."); + } else if (!is_3d && center_words.size() != 2) { + fatal_error("A hexagonal lattice without must have
" + "specified by 2 numbers."); + } + center[0] = stod(center_words[0]); + center[1] = stod(center_words[1]); + if (is_3d) {center[2] = stod(center_words[2]);} + + // Read the lattice pitches. + std::string pitch_str {get_node_value(lat_node, "pitch")}; + std::vector pitch_words {split(pitch_str)}; + if (is_3d && (pitch_words.size() != 2)) { + fatal_error("A hexagonal lattice with must have " + "specified by 2 numbers."); + } else if (!is_3d && (pitch_words.size() != 1)) { + fatal_error("A hexagonal lattice without must have " + "specified by 1 number."); + } + pitch[0] = stod(pitch_words[0]); + if (is_3d) {pitch[1] = stod(pitch_words[1]);} + + // Read the universes and make sure the correct number was specified. + int n_univ = (3*n_rings*n_rings - 3*n_rings + 1) * n_axial; + std::string univ_str {get_node_value(lat_node, "universes")}; + std::vector univ_words {split(univ_str)}; + if (univ_words.size() != n_univ) { + std::stringstream err_msg; + err_msg << "Expected " << n_univ + << " universes for a hexagonal lattice with " << n_rings + << " rings and " << n_axial << " axial levels" << " but " + << univ_words.size() << " were specified."; + fatal_error(err_msg); + } + + // Parse the universes. + // Universes in hexagonal lattices are stored in a manner that represents + // a skewed coordinate system: (x, alpha) rather than (x, y). There is + // no obvious, direct relationship between the order of universes in the + // input and the order that they will be stored in the skewed array so + // the following code walks a set of index values across the skewed array + // in a manner that matches the input order. Note that i_x = 0, i_a = 0 + // corresponds to the center of the hexagonal lattice. + + universes.resize((2*n_rings-1) * (2*n_rings-1) * n_axial, C_NONE); + int input_index = 0; + for (int m = 0; m < n_axial; m++) { + // Initialize lattice indecies. + int i_x = 1; + int i_a = n_rings - 1; + + // Map upper triangular region of hexagonal lattice which is found in the + // first n_rings-1 rows of the input. + for (int k = 0; k < n_rings-1; k++) { + // Walk the index to lower-left neighbor of last row start. + i_x -= 1; + + // Iterate over the input columns. + for (int j = 0; j < k+1; j++) { + int indx = (2*n_rings-1)*(2*n_rings-1) * m + + (2*n_rings-1) * (i_a+n_rings-1) + + (i_x+n_rings-1); + universes[indx] = stoi(univ_words[input_index]); + input_index++; + // Walk the index to the right neighbor (which is not adjacent). + i_x += 2; + i_a -= 1; + } + + // Return the lattice index to the start of the current row. + i_x -= 2 * (k+1); + i_a += (k+1); + } + + // Map the middle square region of the hexagonal lattice which is found in + // the next 2*n_rings-1 rows of the input. + for (int k = 0; k < 2*n_rings-1; k++) { + if ((k % 2) == 0) { + // Walk the index to the lower-left neighbor of the last row start. + i_x -= 1; + } else { + // Walk the index to the lower-right neighbor of the last row start. + i_x += 1; + i_a -= 1; + } + + // Iterate over the input columns. + for (int j = 0; j < n_rings - (k % 2); j++) { + int indx = (2*n_rings-1)*(2*n_rings-1) * m + + (2*n_rings-1) * (i_a+n_rings-1) + + (i_x+n_rings-1); + universes[indx] = stoi(univ_words[input_index]); + input_index++; + // Walk the index to the right neighbor (which is not adjacent). + i_x += 2; + i_a -= 1; + } + + // Return the lattice index to the start of the current row. + i_x -= 2*(n_rings - (k % 2)); + i_a += n_rings - (k % 2); + } + + // Map the lower triangular region of the hexagonal lattice. + for (int k = 0; k < n_rings-1; k++) { + // Walk the index to the lower-right neighbor of the last row start. + i_x += 1; + i_a -= 1; + + // Iterate over the input columns. + for (int j = 0; j < n_rings-k-1; j++) { + int indx = (2*n_rings-1)*(2*n_rings-1) * m + + (2*n_rings-1) * (i_a+n_rings-1) + + (i_x+n_rings-1); + universes[indx] = stoi(univ_words[input_index]); + input_index++; + // Walk the index to the right neighbor (which is not adjacent). + i_x += 2; + i_a -= 1; + } + + // Return lattice index to start of current row. + i_x -= 2*(n_rings - k - 1); + i_a += n_rings - k - 1; + } + } +} + +//============================================================================== + +int32_t& +HexLattice::operator[](const int i_xyz[3]) +{ + int indx = (2*n_rings-1)*(2*n_rings-1) * i_xyz[2] + + (2*n_rings-1) * i_xyz[1] + + i_xyz[0]; + return universes[indx]; +} + +//============================================================================== + +LatticeIter HexLattice::begin() +{return LatticeIter(*this, n_rings-1);} + +ReverseLatticeIter HexLattice::rbegin() +{return ReverseLatticeIter(*this, universes.size()-n_rings);} + +//============================================================================== + +bool +HexLattice::are_valid_indices(const int i_xyz[3]) const +{ + return ((i_xyz[0] >= 0) && (i_xyz[1] >= 0) && (i_xyz[2] >= 0) + && (i_xyz[0] < 2*n_rings-1) && (i_xyz[1] < 2*n_rings-1) + && (i_xyz[0] + i_xyz[1] > n_rings-2) + && (i_xyz[0] + i_xyz[1] < 3*n_rings-2) + && (i_xyz[2] < n_axial)); +} + +//============================================================================== + +std::pair> +HexLattice::distance(const double xyz[3], const double uvw[3], + const int i_xyz[3]) const +{ + // Compute the direction on the hexagonal basis. + double beta_dir = uvw[0] * std::sqrt(3.0) / 2.0 + uvw[1] / 2.0; + double gamma_dir = uvw[0] * std::sqrt(3.0) / 2.0 - uvw[1] / 2.0; + + // Note that hexagonal lattice distance calculations are performed + // using the particle's coordinates relative to the neighbor lattice + // cells, not relative to the particle's current cell. This is done + // because there is significant disagreement between neighboring cells + // on where the lattice boundary is due to finite precision issues. + + // Upper-right and lower-left sides. + double d {INFTY}; + std::array lattice_trans; + double edge = -copysign(0.5*pitch[0], beta_dir); // Oncoming edge + std::array xyz_t; + if (beta_dir > 0) { + const int i_xyz_t[3] {i_xyz[0]+1, i_xyz[1], i_xyz[2]}; + xyz_t = get_local_xyz(xyz, i_xyz_t); + } else { + const int i_xyz_t[3] {i_xyz[0]-1, i_xyz[1], i_xyz[2]}; + xyz_t = get_local_xyz(xyz, i_xyz_t); + } + double beta = xyz_t[0] * std::sqrt(3.0) / 2.0 + xyz_t[1] / 2.0; + if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) { + d = (edge - beta) / beta_dir; + if (beta_dir > 0) { + lattice_trans = {1, 0, 0}; + } else { + lattice_trans = {-1, 0, 0}; + } + } + + // Lower-right and upper-left sides. + edge = -copysign(0.5*pitch[0], gamma_dir); + if (gamma_dir > 0) { + const int i_xyz_t[3] {i_xyz[0]+1, i_xyz[1]-1, i_xyz[2]}; + xyz_t = get_local_xyz(xyz, i_xyz_t); + } else { + const int i_xyz_t[3] {i_xyz[0]-1, i_xyz[1]+1, i_xyz[2]}; + xyz_t = get_local_xyz(xyz, i_xyz_t); + } + double gamma = xyz_t[0] * std::sqrt(3.0) / 2.0 - xyz_t[1] / 2.0; + if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) { + double this_d = (edge - gamma) / gamma_dir; + if (this_d < d) { + if (gamma_dir > 0) { + lattice_trans = {1, -1, 0}; + } else { + lattice_trans = {-1, 1, 0}; + } + d = this_d; + } + } + + // Upper and lower sides. + edge = -copysign(0.5*pitch[0], uvw[1]); + if (uvw[1] > 0) { + const int i_xyz_t[3] {i_xyz[0], i_xyz[1]+1, i_xyz[2]}; + xyz_t = get_local_xyz(xyz, i_xyz_t); + } else { + const int i_xyz_t[3] {i_xyz[0], i_xyz[1]-1, i_xyz[2]}; + xyz_t = get_local_xyz(xyz, i_xyz_t); + } + if ((std::abs(xyz_t[1] - edge) > FP_PRECISION) && uvw[1] != 0) { + double this_d = (edge - xyz_t[1]) / uvw[1]; + if (this_d < d) { + if (uvw[1] > 0) { + lattice_trans = {0, 1, 0}; + } else { + lattice_trans = {0, -1, 0}; + } + d = this_d; + } + } + + // Top and bottom sides + if (is_3d) { + double z {xyz[2]}; + double w {uvw[2]}; + double z0 {copysign(0.5 * pitch[1], w)}; + if ((std::abs(z - z0) > FP_PRECISION) && w != 0) { + double this_d = (z0 - z) / w; + if (this_d < d) { + d = this_d; + if (w > 0) { + lattice_trans = {0, 0, 1}; + } else { + lattice_trans = {0, 0, -1}; + } + d = this_d; + } + } + } + + return {d, lattice_trans}; +} + +//============================================================================== + +std::array +HexLattice::get_indices(const double xyz[3]) const +{ + // Offset the xyz by the lattice center. + double xyz_o[3] {xyz[0] - center[0], xyz[1] - center[1], xyz[2]}; + if (is_3d) {xyz_o[2] -= center[2];} + + // Index the z direction. + std::array out; + if (is_3d) { + out[2] = static_cast(std::ceil(xyz_o[2] / pitch[1] + 0.5 * n_axial))-1; + } else { + out[2] = 0; + } + + // Convert coordinates into skewed bases. The (x, alpha) basis is used to + // find the index of the global coordinates to within 4 cells. + double alpha = xyz_o[1] - xyz_o[0] / std::sqrt(3.0); + out[0] = static_cast(std::floor(xyz_o[0] + / (0.5*std::sqrt(3.0) * pitch[0]))); + out[1] = static_cast(std::floor(alpha / pitch[0])); + + // Add offset to indices (the center cell is (i_x, i_alpha) = (0, 0) but + // the array is offset so that the indices never go below 0). + out[0] += n_rings-1; + out[1] += n_rings-1; + + // Calculate the (squared) distance between the particle and the centers of + // the four possible cells. Regular hexagonal tiles form a Voronoi + // tessellation so the xyz should be in the hexagonal cell that it is closest + // to the center of. This method is used over a method that uses the + // remainders of the floor divisions above because it provides better finite + // precision performance. Squared distances are used becasue they are more + // computationally efficient than normal distances. + int k {1}; + int k_min {1}; + double d_min {INFTY}; + for (int i = 0; i < 2; i++) { + for (int j = 0; j < 2; j++) { + int i_xyz[3] {out[0] + j, out[1] + i, 0}; + std::array xyz_t = get_local_xyz(xyz, i_xyz); + double d = xyz_t[0]*xyz_t[0] + xyz_t[1]*xyz_t[1]; + if (d < d_min) { + d_min = d; + k_min = k; + } + k++; + } + } + + // Select the minimum squared distance which corresponds to the cell the + // coordinates are in. + if (k_min == 2) { + out[0] += 1; + } else if (k_min == 3) { + out[1] += 1; + } else if (k_min == 4) { + out[0] += 1; + out[1] += 1; + } + + return out; +} + +//============================================================================== + +std::array +HexLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const +{ + std::array local_xyz; + + // x_l = x_g - (center + pitch_x*cos(30)*index_x) + local_xyz[0] = global_xyz[0] - (center[0] + + std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings + 1) * pitch[0]); + // y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y) + local_xyz[1] = global_xyz[1] - (center[1] + + (i_xyz[1] - n_rings + 1) * pitch[0] + + (i_xyz[0] - n_rings + 1) * pitch[0] / 2.0); + if (is_3d) { + local_xyz[2] = global_xyz[2] - center[2] + + (0.5 * n_axial - i_xyz[2] - 0.5) * pitch[1]; + } else { + local_xyz[2] = global_xyz[2]; + } + + return local_xyz; +} + +//============================================================================== + +bool +HexLattice::is_valid_index(int indx) const +{ + int nx {2*n_rings - 1}; + int ny {2*n_rings - 1}; + int nz {n_axial}; + int iz = indx / (nx * ny); + int iy = (indx - nx*ny*iz) / nx; + int ix = indx - nx*ny*iz - nx*iy; + int i_xyz[3] {ix, iy, iz}; + return are_valid_indices(i_xyz); +} + +//============================================================================== + +int32_t& +HexLattice::offset(int map, const int i_xyz[3]) +{ + int nx {2*n_rings - 1}; + int ny {2*n_rings - 1}; + int nz {n_axial}; + return offsets[nx*ny*nz*map + nx*ny*i_xyz[2] + nx*i_xyz[1] + i_xyz[0]]; +} + +//============================================================================== + +std::string +HexLattice::index_to_string(int indx) const +{ + int nx {2*n_rings - 1}; + int ny {2*n_rings - 1}; + int iz {indx / (nx * ny)}; + int iy {(indx - nx * ny * iz) / nx}; + int ix {indx - nx * ny * iz - nx * iy}; + std::string out {std::to_string(ix - n_rings + 1)}; + out += ','; + out += std::to_string(iy - n_rings + 1); + if (is_3d) { + out += ','; + out += std::to_string(iz); + } + return out; +} + +//============================================================================== + +void +HexLattice::to_hdf5_inner(hid_t lat_group) const +{ + // Write basic lattice information. + write_string(lat_group, "type", "hexagonal", false); + write_int(lat_group, 0, nullptr, "n_rings", &n_rings, false); + write_int(lat_group, 0, nullptr, "n_axial", &n_axial, false); + if (is_3d) { + write_double(lat_group, "pitch", pitch, false); + write_double(lat_group, "center", center, false); + } else { + std::array pitch_short {{pitch[0]}}; + write_double(lat_group, "pitch", pitch_short, false); + std::array center_short {{center[0], center[1]}}; + write_double(lat_group, "center", center_short, false); + } + + // Write the universe ids. + hsize_t nx {static_cast(2*n_rings - 1)}; + hsize_t ny {static_cast(2*n_rings - 1)}; + hsize_t nz {static_cast(n_axial)}; + int out[nx*ny*nz]; + + for (int m = 0; m < nz; m++) { + for (int k = 0; k < ny; k++) { + for (int j = 0; j < nx; j++) { + int indx = nx*ny*m + nx*k + j; + if (j + k < n_rings - 1) { + // This array position is never used; put a -1 to indicate this. + out[indx] = -1; + } else if (j + k > 3*n_rings - 3) { + // This array position is never used; put a -1 to indicate this. + out[indx] = -1; + } else { + out[indx] = global_universes[universes[indx]]->id; + } + } + } + } + + hsize_t dims[3] {nz, ny, nx}; + write_int(lat_group, 3, dims, "universes", out, false); +} + +//============================================================================== +// Non-method functions +//============================================================================== + +extern "C" void +read_lattices(pugi::xml_node *node) +{ + for (pugi::xml_node lat_node : node->children("lattice")) { + lattices_c.push_back(new RectLattice(lat_node)); + } + for (pugi::xml_node lat_node : node->children("hex_lattice")) { + lattices_c.push_back(new HexLattice(lat_node)); + } + + // Fill the lattice map. + for (int i_lat = 0; i_lat < lattices_c.size(); i_lat++) { + int id = lattices_c[i_lat]->id; + auto in_map = lattice_map.find(id); + if (in_map == lattice_map.end()) { + lattice_map[id] = i_lat; + } else { + std::stringstream err_msg; + err_msg << "Two or more lattices use the same unique ID: " << id; + fatal_error(err_msg); + } + } +} + +//============================================================================== +// Fortran compatibility functions +//============================================================================== + +extern "C" { + Lattice* lattice_pointer(int lat_ind) {return lattices_c[lat_ind];} + + int32_t lattice_id(Lattice *lat) {return lat->id;} + + bool lattice_are_valid_indices(Lattice *lat, const int i_xyz[3]) + {return lat->are_valid_indices(i_xyz);} + + void lattice_distance(Lattice *lat, const double xyz[3], const double uvw[3], + const int i_xyz[3], double *d, int lattice_trans[3]) + { + std::pair> ld {lat->distance(xyz, uvw, i_xyz)}; + *d = ld.first; + lattice_trans[0] = ld.second[0]; + lattice_trans[1] = ld.second[1]; + lattice_trans[2] = ld.second[2]; + } + + void lattice_get_indices(Lattice *lat, const double xyz[3], int i_xyz[3]) + { + std::array inds = lat->get_indices(xyz); + i_xyz[0] = inds[0]; + i_xyz[1] = inds[1]; + i_xyz[2] = inds[2]; + } + + void lattice_get_local_xyz(Lattice *lat, const double global_xyz[3], + const int i_xyz[3], double local_xyz[3]) + { + std::array xyz = lat->get_local_xyz(global_xyz, i_xyz); + local_xyz[0] = xyz[0]; + local_xyz[1] = xyz[1]; + local_xyz[2] = xyz[2]; + } + + int32_t lattice_offset(Lattice *lat, int map, const int i_xyz[3]) + {return lat->offset(map, i_xyz);} + + int32_t lattice_outer(Lattice *lat) {return lat->outer;} + + void lattice_to_hdf5(Lattice *lat, hid_t group) {lat->to_hdf5(group);} + + int32_t lattice_universe(Lattice *lat, const int i_xyz[3]) + {return (*lat)[i_xyz];} +} + +} // namespace openmc diff --git a/src/lattice.h b/src/lattice.h new file mode 100644 index 0000000000..19749ce4f0 --- /dev/null +++ b/src/lattice.h @@ -0,0 +1,259 @@ +#ifndef LATTICE_H +#define LATTICE_H + +#include +#include +#include +#include +#include + +#include "constants.h" +#include "hdf5.h" +#include "pugixml/pugixml.hpp" + + +namespace openmc { + +//============================================================================== +// Module constants +//============================================================================== + +constexpr int32_t NO_OUTER_UNIVERSE{-1}; + +//============================================================================== +// Global variables +//============================================================================== + +class Lattice; +extern std::vector lattices_c; + +extern std::unordered_map lattice_map; + +//============================================================================== +//! \class Lattice +//! \brief Abstract type for ordered array of universes. +//============================================================================== + +class LatticeIter; +class ReverseLatticeIter; + +class Lattice +{ +public: + int32_t id; //!< Universe ID number + std::string name; //!< User-defined name + std::vector universes; //!< Universes filling each lattice tile + int32_t outer {NO_OUTER_UNIVERSE}; //!< Universe tiled outside the lattice + std::vector offsets; //!< Distribcell offset table + + explicit Lattice(pugi::xml_node lat_node); + + virtual ~Lattice() {} + + virtual int32_t& operator[](const int i_xyz[3]) = 0; + + virtual LatticeIter begin(); + LatticeIter end(); + + virtual ReverseLatticeIter rbegin(); + ReverseLatticeIter rend(); + + //! Convert internal universe values from IDs to indices using universe_map. + void adjust_indices(); + + //! Allocate offset table for distribcell. + void allocate_offset_table(int n_maps) + {offsets.resize(n_maps * universes.size(), C_NONE);} + + //! Populate the distribcell offset tables. + int32_t fill_offset_table(int32_t offset, int32_t target_univ_id, int map); + + //! \brief Check lattice indices. + //! @param i_xyz[3] The indices for a lattice tile. + //! @return true if the given indices fit within the lattice bounds. False + //! otherwise. + virtual bool are_valid_indices(const int i_xyz[3]) const = 0; + + //! \brief Find the next lattice surface crossing + //! @param xyz[3] A 3D Cartesian coordinate. + //! @param uvw[3] A 3D Cartesian direction. + //! @param i_xyz[3] The indices for a lattice tile. + //! @return The distance to the next crossing and an array indicating how the + //! lattice indices would change after crossing that boundary. + virtual std::pair> + distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const + = 0; + + //! \brief Find the lattice tile indices for a given point. + //! @param xyz[3] A 3D Cartesian coordinate. + //! @return An array containing the indices of a lattice tile. + virtual std::array get_indices(const double xyz[3]) const = 0; + + //! \brief Get coordinates local to a lattice tile. + //! @param global_xyz[3] A 3D Cartesian coordinate. + //! @param i_xyz[3] The indices for a lattice tile. + //! @return Local 3D Cartesian coordinates. + virtual std::array + get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const = 0; + + //! \brief Check flattened lattice index. + //! @param indx The index for a lattice tile. + //! @return true if the given index fit within the lattice bounds. False + //! otherwise. + virtual bool is_valid_index(int indx) const + {return (indx >= 0) && (indx < universes.size());} + + //! \brief Get the distribcell offset for a lattice tile. + //! @param The map index for the target cell. + //! @param i_xyz[3] The indices for a lattice tile. + //! @return Distribcell offset i.e. the largest instance number for the target + //! cell found in the geometry tree under this lattice tile. + virtual int32_t& offset(int map, const int i_xyz[3]) = 0; + + //! \brief Convert an array index to a useful human-readable string. + //! @param indx The index for a lattice tile. + //! @return A string representing the lattice tile. + virtual std::string index_to_string(int indx) const = 0; + + //! \brief Write lattice information to an HDF5 group. + //! @param group_id An HDF5 group id. + void to_hdf5(hid_t group_id) const; + +protected: + bool is_3d; //!< Has divisions along the z-axis? + + virtual void to_hdf5_inner(hid_t group_id) const = 0; +}; + +//============================================================================== +//! An iterator over lattice universes. +//============================================================================== + +class LatticeIter +{ +public: + int indx; //!< An index to a Lattice universes or offsets array. + + LatticeIter(Lattice &lat_, int indx_) + : lat(lat_), + indx(indx_) + {} + + bool operator==(const LatticeIter &rhs) {return (indx == rhs.indx);} + + bool operator!=(const LatticeIter &rhs) {return !(*this == rhs);} + + int32_t& operator*() {return lat.universes[indx];} + + LatticeIter& operator++() + { + while (indx < lat.universes.size()) { + ++indx; + if (lat.is_valid_index(indx)) return *this; + } + indx = lat.universes.size(); + return *this; + } + +protected: + Lattice ⪫ +}; + +//============================================================================== +//! A reverse iterator over lattice universes. +//============================================================================== + +class ReverseLatticeIter : public LatticeIter +{ +public: + ReverseLatticeIter(Lattice &lat_, int indx_) + : LatticeIter {lat_, indx_} + {} + + ReverseLatticeIter& operator++() + { + while (indx > -1) { + --indx; + if (lat.is_valid_index(indx)) return *this; + } + indx = -1; + return *this; + } +}; + +//============================================================================== + +class RectLattice : public Lattice +{ +public: + explicit RectLattice(pugi::xml_node lat_node); + + int32_t& operator[](const int i_xyz[3]); + + bool are_valid_indices(const int i_xyz[3]) const; + + std::pair> + distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const; + + std::array get_indices(const double xyz[3]) const; + + std::array + get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const; + + int32_t& offset(int map, const int i_xyz[3]); + + std::string index_to_string(int indx) const; + + void to_hdf5_inner(hid_t group_id) const; + +private: + std::array n_cells; //!< Number of cells along each axis + std::array lower_left; //!< Global lower-left corner of the lattice + std::array pitch; //!< Lattice tile width along each axis + + // Convenience aliases + int &nx {n_cells[0]}; + int &ny {n_cells[1]}; + int &nz {n_cells[2]}; +}; + +//============================================================================== + +class HexLattice : public Lattice +{ +public: + explicit HexLattice(pugi::xml_node lat_node); + + int32_t& operator[](const int i_xyz[3]); + + LatticeIter begin(); + + ReverseLatticeIter rbegin(); + + bool are_valid_indices(const int i_xyz[3]) const; + + std::pair> + distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const; + + std::array get_indices(const double xyz[3]) const; + + std::array + get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const; + + bool is_valid_index(int indx) const; + + int32_t& offset(int map, const int i_xyz[3]); + + std::string index_to_string(int indx) const; + + void to_hdf5_inner(hid_t group_id) const; + +private: + int n_rings; //!< Number of radial tile positions + int n_axial; //!< Number of axial tile positions + std::array center; //!< Global center of lattice + std::array pitch; //!< Lattice tile width and height +}; + +} // namespace openmc +#endif // LATTICE_H diff --git a/src/mgxs.cpp b/src/mgxs.cpp index 3a079a9325..ab2970b2f3 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -11,7 +11,8 @@ std::vector macro_xs; // Mgxs base-class methods //============================================================================== -void Mgxs::init(const std::string& in_name, const double in_awr, +void +Mgxs::init(const std::string& in_name, const double in_awr, const double_1dvec& in_kTs, const bool in_fissionable, const int in_scatter_format, const int in_num_groups, const int in_num_delayed_groups, const double_1dvec& in_polar, @@ -41,8 +42,10 @@ void Mgxs::init(const std::string& in_name, const double in_awr, } } +//============================================================================== -void Mgxs::_metadata_from_hdf5(const hid_t xs_id, const int in_num_groups, +void +Mgxs::_metadata_from_hdf5(const hid_t xs_id, const int in_num_groups, const int in_num_delayed_groups, double_1dvec& temperature, int& method, const double tolerance, int_1dvec& temps_to_read, int& order_dim, bool& is_isotropic, const int n_threads) @@ -249,8 +252,10 @@ void Mgxs::_metadata_from_hdf5(const hid_t xs_id, const int in_num_groups, n_threads); } +//============================================================================== -void Mgxs::from_hdf5(hid_t xs_id, const int energy_groups, +void +Mgxs::from_hdf5(hid_t xs_id, const int energy_groups, const int delayed_groups, double_1dvec& temperature, int& method, const double tolerance, const int max_order, const bool legendre_to_tabular, const int legendre_to_tabular_points, @@ -288,8 +293,10 @@ void Mgxs::from_hdf5(hid_t xs_id, const int energy_groups, scatter_format = final_scatter_format; } +//============================================================================== -void Mgxs::build_macro(const std::string& in_name, double_1dvec& mat_kTs, +void +Mgxs::build_macro(const std::string& in_name, double_1dvec& mat_kTs, std::vector& micros, double_1dvec& atom_densities, int& method, const double tolerance, const int n_threads) { @@ -384,9 +391,11 @@ void Mgxs::build_macro(const std::string& in_name, double_1dvec& mat_kTs, } // end temperature (t) loop } +//============================================================================== -void Mgxs::combine(std::vector& micros, double_1dvec& scalars, - int_1dvec& micro_ts, int this_t) +void +Mgxs::combine(std::vector& micros, double_1dvec& scalars, + int_1dvec& micro_ts, int this_t) { // Build the vector of pointers to the xs objects within micros std::vector those_xs(micros.size()); @@ -400,9 +409,11 @@ void Mgxs::combine(std::vector& micros, double_1dvec& scalars, xs[this_t].combine(those_xs, scalars); } +//============================================================================== -double Mgxs::get_xs(const int tid, const int xstype, const int gin, - int* gout, double* mu, int* dg) +double +Mgxs::get_xs(const int tid, const int xstype, const int gin, int* gout, + double* mu, int* dg) { // This method assumes that the temperature and angle indices are set double val; @@ -524,9 +535,10 @@ double Mgxs::get_xs(const int tid, const int xstype, const int gin, return val; } +//============================================================================== -void Mgxs::sample_fission_energy(const int tid, const int gin, int& dg, - int& gout) +void +Mgxs::sample_fission_energy(const int tid, const int gin, int& dg, int& gout) { // This method assumes that the temperature and angle indices are set double nu_fission = xs[cache[tid].t].nu_fission[cache[tid].p][cache[tid].a][gin]; @@ -581,8 +593,10 @@ void Mgxs::sample_fission_energy(const int tid, const int gin, int& dg, } } +//============================================================================== -void Mgxs::sample_scatter(const int tid, const int gin, int& gout, double& mu, +void +Mgxs::sample_scatter(const int tid, const int gin, int& gout, double& mu, double& wgt) { // This method assumes that the temperature and angle indices are set @@ -590,8 +604,10 @@ void Mgxs::sample_scatter(const int tid, const int gin, int& gout, double& mu, xs[cache[tid].t].scatter[cache[tid].p][cache[tid].a]->sample(gin, gout, mu, wgt); } +//============================================================================== -void Mgxs::calculate_xs(const int tid, const int gin, const double sqrtkT, +void +Mgxs::calculate_xs(const int tid, const int gin, const double sqrtkT, const double uvw[3], double& total_xs, double& abs_xs, double& nu_fiss_xs) { // Set our indices @@ -607,8 +623,10 @@ void Mgxs::calculate_xs(const int tid, const int gin, const double sqrtkT, } } +//============================================================================== -bool Mgxs::equiv(const Mgxs& that) +bool +Mgxs::equiv(const Mgxs& that) { bool match = false; @@ -624,8 +642,10 @@ bool Mgxs::equiv(const Mgxs& that) return match; } +//============================================================================== -void Mgxs::set_temperature_index(const int tid, const double sqrtkT) +void +Mgxs::set_temperature_index(const int tid, const double sqrtkT) { // See if we need to find the new index if (sqrtkT != cache[tid].sqrtkT) { @@ -644,8 +664,10 @@ void Mgxs::set_temperature_index(const int tid, const double sqrtkT) } } +//============================================================================== -void Mgxs::set_angle_index(const int tid, const double uvw[3]) +void +Mgxs::set_angle_index(const int tid, const double uvw[3]) { // See if we need to find the new index if ((uvw[0] != cache[tid].uvw[0]) || (uvw[1] != cache[tid].uvw[1]) || diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp index 6c9dc9f42e..29202e49d5 100644 --- a/src/mgxs_interface.cpp +++ b/src/mgxs_interface.cpp @@ -6,7 +6,8 @@ namespace openmc { // Mgxs data loading interface methods //============================================================================== -void add_mgxs_c(hid_t file_id, char* name, const int energy_groups, +void +add_mgxs_c(hid_t file_id, char* name, const int energy_groups, const int delayed_groups, const int n_temps, double temps[], int& method, const double tolerance, const int max_order, const bool legendre_to_tabular, const int legendre_to_tabular_points, @@ -38,8 +39,10 @@ void add_mgxs_c(hid_t file_id, char* name, const int energy_groups, nuclides_MG.push_back(mg); } +//============================================================================== -bool query_fissionable_c(const int n_nuclides, const int i_nuclides[]) +bool +query_fissionable_c(const int n_nuclides, const int i_nuclides[]) { bool result = false; for (int n = 0; n < n_nuclides; n++) { @@ -48,8 +51,10 @@ bool query_fissionable_c(const int n_nuclides, const int i_nuclides[]) return result; } +//============================================================================== -void create_macro_xs_c(char* mat_name, const int n_nuclides, +void +create_macro_xs_c(char* mat_name, const int n_nuclides, const int i_nuclides[], const int n_temps, const double temps[], const double atom_densities[], int& method, const double tolerance, const int n_threads) @@ -81,7 +86,8 @@ void create_macro_xs_c(char* mat_name, const int n_nuclides, // Mgxs tracking/transport/tallying interface methods //============================================================================== -void calculate_xs_c(const int i_mat, const int tid, const int gin, +void +calculate_xs_c(const int i_mat, const int tid, const int gin, const double sqrtkT, const double uvw[3], double& total_xs, double& abs_xs, double& nu_fiss_xs) { @@ -89,8 +95,10 @@ void calculate_xs_c(const int i_mat, const int tid, const int gin, nu_fiss_xs); } +//============================================================================== -void sample_scatter_c(const int i_mat, const int tid, const int gin, int& gout, +void +sample_scatter_c(const int i_mat, const int tid, const int gin, int& gout, double& mu, double& wgt, double uvw[3]) { int gout_c = gout - 1; @@ -103,8 +111,10 @@ void sample_scatter_c(const int i_mat, const int tid, const int gin, int& gout, rotate_angle_c(uvw, mu, nullptr); } +//============================================================================== -void sample_fission_energy_c(const int i_mat, const int tid, const int gin, +void +sample_fission_energy_c(const int i_mat, const int tid, const int gin, int& dg, int& gout) { int dg_c = 0; @@ -116,8 +126,10 @@ void sample_fission_energy_c(const int i_mat, const int tid, const int gin, gout = gout_c + 1; } +//============================================================================== -double get_nuclide_xs_c(const int index, const int tid, const int xstype, +double +get_nuclide_xs_c(const int index, const int tid, const int xstype, const int gin, int* gout, double* mu, int* dg) { int gout_c; @@ -140,8 +152,10 @@ double get_nuclide_xs_c(const int index, const int tid, const int xstype, dg_c_p); } +//============================================================================== -double get_macro_xs_c(const int index, const int tid, const int xstype, +double +get_macro_xs_c(const int index, const int tid, const int xstype, const int gin, int* gout, double* mu, int* dg) { int gout_c; @@ -164,36 +178,42 @@ double get_macro_xs_c(const int index, const int tid, const int xstype, dg_c_p); } +//============================================================================== -void set_nuclide_angle_index_c(const int index, const int tid, +void +set_nuclide_angle_index_c(const int index, const int tid, const double uvw[3]) { // Update the values nuclides_MG[index - 1].set_angle_index(tid, uvw); } +//============================================================================== -void set_macro_angle_index_c(const int index, const int tid, +void +set_macro_angle_index_c(const int index, const int tid, const double uvw[3]) { // Update the values macro_xs[index - 1].set_angle_index(tid, uvw); } +//============================================================================== -void set_nuclide_temperature_index_c(const int index, const int tid, +void +set_nuclide_temperature_index_c(const int index, const int tid, const double sqrtkT) { // Update the values nuclides_MG[index - 1].set_temperature_index(tid, sqrtkT); } - //============================================================================== -// Mgxs general methods +// General Mgxs methods //============================================================================== -void get_name_c(const int index, int name_len, char* name) +void +get_name_c(const int index, int name_len, char* name) { // First blank out our input string std::string str(name_len, ' '); @@ -207,8 +227,10 @@ void get_name_c(const int index, int name_len, char* name) name[std::strlen(name)] = ' '; } +//============================================================================== -double get_awr_c(const int index) +double +get_awr_c(const int index) { return nuclides_MG[index - 1].awr; } diff --git a/src/mgxs_interface.h b/src/mgxs_interface.h index 4e10682669..09993e9811 100644 --- a/src/mgxs_interface.h +++ b/src/mgxs_interface.h @@ -9,51 +9,79 @@ namespace openmc { +//============================================================================== +// Global variables +//============================================================================== + extern std::vector nuclides_MG; extern std::vector macro_xs; +//============================================================================== +// Mgxs data loading interface methods +//============================================================================== -extern "C" void add_mgxs_c(hid_t file_id, char* name, const int energy_groups, +extern "C" void +add_mgxs_c(hid_t file_id, char* name, const int energy_groups, const int delayed_groups, const int n_temps, double temps[], int& method, const double tolerance, const int max_order, const bool legendre_to_tabular, const int legendre_to_tabular_points, const int n_threads); -extern "C" bool query_fissionable_c(const int n_nuclides, const int i_nuclides[]); +extern "C" bool +query_fissionable_c(const int n_nuclides, const int i_nuclides[]); -extern "C" void create_macro_xs_c(char* mat_name, const int n_nuclides, +extern "C" void +create_macro_xs_c(char* mat_name, const int n_nuclides, const int i_nuclides[], const int n_temps, const double temps[], const double atom_densities[], int& method, const double tolerance, const int n_threads); -extern "C" void calculate_xs_c(const int i_mat, const int tid, const int gin, +//============================================================================== +// Mgxs tracking/transport/tallying interface methods +//============================================================================== + +extern "C" void +calculate_xs_c(const int i_mat, const int tid, const int gin, const double sqrtkT, const double uvw[3], double& total_xs, double& abs_xs, double& nu_fiss_xs); -extern "C" void sample_scatter_c(const int i_mat, const int tid, const int gin, +extern "C" void +sample_scatter_c(const int i_mat, const int tid, const int gin, int& gout, double& mu, double& wgt, double uvw[3]); -extern "C" void sample_fission_energy_c(const int i_mat, const int tid, +extern "C" void +sample_fission_energy_c(const int i_mat, const int tid, const int gin, int& dg, int& gout); -extern "C" void get_name_c(const int index, int name_len, char* name); - -extern "C" double get_awr_c(const int index); - -extern "C" double get_nuclide_xs_c(const int index, const int tid, +extern "C" double +get_nuclide_xs_c(const int index, const int tid, const int xstype, const int gin, int* gout, double* mu, int* dg); -extern "C" double get_macro_xs_c(const int index, const int tid, +extern "C" double +get_macro_xs_c(const int index, const int tid, const int xstype, const int gin, int* gout, double* mu, int* dg); -extern "C" void set_nuclide_angle_index_c(const int index, const int tid, +extern "C" void +set_nuclide_angle_index_c(const int index, const int tid, const double uvw[3]); -extern "C" void set_macro_angle_index_c(const int index, const int tid, +extern "C" void +set_macro_angle_index_c(const int index, const int tid, const double uvw[3]); -extern "C" void set_nuclide_temperature_index_c(const int index, const int tid, +extern "C" void +set_nuclide_temperature_index_c(const int index, const int tid, const double sqrtkT); +//============================================================================== +// General Mgxs methods +//============================================================================== + +extern "C" void +get_name_c(const int index, int name_len, char* name); + +extern "C" double +get_awr_c(const int index); + } // namespace openmc #endif // MGXS_INTERFACE_H \ No newline at end of file diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 94aa1928ad..46f5b913ed 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -10,7 +10,6 @@ module nuclide_header use endf_header, only: Function1D, Polynomial, Tabulated1D use error use hdf5_interface - use list_header, only: ListInt use math, only: faddeeva, w_derivative, & broaden_wmp_polynomials use multipole_header, only: FORM_RM, FORM_MLBW, MP_EA, RM_RT, RM_RA, & diff --git a/src/output.F90 b/src/output.F90 index 9c5d9cb451..6197cdc11c 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -234,7 +234,7 @@ contains ! Print cell for this level if (p % coord(i) % cell /= NONE) then c => cells(p % coord(i) % cell) - write(ou,*) ' Cell = ' // trim(to_str(c % id)) + write(ou,*) ' Cell = ' // trim(to_str(c % id())) end if ! Print universe for this level @@ -246,7 +246,7 @@ contains ! Print information on lattice if (p % coord(i) % lattice /= NONE) then l => lattices(p % coord(i) % lattice) % obj - write(ou,*) ' Lattice = ' // trim(to_str(l % id)) + write(ou,*) ' Lattice = ' // trim(to_str(l % id())) write(ou,*) ' Lattice position = (' // trim(to_str(& p % coord(i) % lattice_x)) // ',' // trim(to_str(& p % coord(i) % lattice_y)) // ')' @@ -258,7 +258,7 @@ contains end do ! Print surface - if (p % surface /= NONE) then + if (p % surface /= ERROR_INT) then write(ou,*) ' Surface = ' // to_str(sign(surfaces(i)%id(), p % surface)) end if @@ -633,7 +633,7 @@ contains write(ou,100) 'Cell ID','No. Overlap Checks' do i = 1, n_cells - write(ou,101) cells(i) % id, overlap_check_cnt(i) + write(ou,101) cells(i) % id(), overlap_check_cnt(i) if (overlap_check_cnt(i) < 10) num_sparse = num_sparse + 1 end do write(ou,*) @@ -643,7 +643,7 @@ contains do i = 1, n_cells if (overlap_check_cnt(i) < 10) then j = j + 1 - write(ou,'(1X,A8)', advance='no') trim(to_str(cells(i) % id)) + write(ou,'(1X,A8)', advance='no') trim(to_str(cells(i) % id())) if (modulo(j,8) == 0) write(ou,*) end if end do diff --git a/src/plot.F90 b/src/plot.F90 index f17ac33511..789ad0ee9b 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -85,7 +85,7 @@ contains if (pl % color_by == PLOT_COLOR_MATS) then ! Assign color based on material associate (c => cells(p % coord(j) % cell)) - if (c % type == FILL_UNIVERSE) then + if (c % type() == FILL_UNIVERSE) then ! If we stopped on a middle universe level, treat as if not found rgb = pl % not_found % rgb id = -1 @@ -101,7 +101,7 @@ contains else if (pl % color_by == PLOT_COLOR_CELLS) then ! Assign color based on cell rgb = pl % colors(p % coord(j) % cell) % rgb - id = cells(p % coord(j) % cell) % id + id = cells(p % coord(j) % cell) % id() else rgb = 0 id = -1 diff --git a/src/scattdata.cpp b/src/scattdata.cpp index 1e33c07c17..ad6a1db697 100644 --- a/src/scattdata.cpp +++ b/src/scattdata.cpp @@ -6,8 +6,9 @@ namespace openmc { // ScattData base-class methods //============================================================================== -void ScattData::generic_init(int order, int_1dvec& in_gmin, - int_1dvec& in_gmax, double_2dvec& in_energy, double_2dvec& in_mult) +void +ScattData::generic_init(int order, int_1dvec& in_gmin, int_1dvec& in_gmax, + double_2dvec& in_energy, double_2dvec& in_mult) { int groups = in_energy.size(); @@ -38,8 +39,10 @@ void ScattData::generic_init(int order, int_1dvec& in_gmin, } } +//============================================================================== -void ScattData::sample_energy(int gin, int& gout, int& i_gout) +void +ScattData::sample_energy(int gin, int& gout, int& i_gout) { // Sample the outgoing group double xi = prn(); @@ -54,8 +57,10 @@ void ScattData::sample_energy(int gin, int& gout, int& i_gout) } } +//============================================================================== -double ScattData::get_xs(const int xstype, int gin, int* gout, double* mu) +double +ScattData::get_xs(const int xstype, int gin, int* gout, double* mu) { // Set the outgoing group offset index as needed int i_gout = 0; @@ -109,13 +114,13 @@ double ScattData::get_xs(const int xstype, int gin, int* gout, double* mu) return val; } - //============================================================================== // ScattDataLegendre methods //============================================================================== -void ScattDataLegendre::init(int_1dvec& in_gmin, int_1dvec& in_gmax, - double_2dvec& in_mult, double_3dvec& coeffs) +void +ScattDataLegendre::init(int_1dvec& in_gmin, int_1dvec& in_gmax, + double_2dvec& in_mult, double_3dvec& coeffs) { int groups = coeffs.size(); int order = coeffs[0][0].size(); @@ -169,8 +174,10 @@ void ScattDataLegendre::init(int_1dvec& in_gmin, int_1dvec& in_gmax, update_max_val(); } +//============================================================================== -void ScattDataLegendre::update_max_val() +void +ScattDataLegendre::update_max_val() { int groups = max_val.size(); // Step through the polynomial with fixed number of points to identify the @@ -204,8 +211,10 @@ void ScattDataLegendre::update_max_val() } } +//============================================================================== -double ScattDataLegendre::calc_f(int gin, int gout, double mu) +double +ScattDataLegendre::calc_f(int gin, int gout, double mu) { double f; if ((gout < gmin[gin]) || (gout > gmax[gin])) { @@ -218,8 +227,10 @@ double ScattDataLegendre::calc_f(int gin, int gout, double mu) return f; } +//============================================================================== -void ScattDataLegendre::sample(int gin, int& gout, double& mu, double& wgt) +void +ScattDataLegendre::sample(int gin, int& gout, double& mu, double& wgt) { // Sample the outgoing energy using the base-class method int i_gout; @@ -247,9 +258,11 @@ void ScattDataLegendre::sample(int gin, int& gout, double& mu, double& wgt) wgt *= mult[gin][i_gout]; } +//============================================================================== -void ScattDataLegendre::combine(std::vector& those_scatts, - double_1dvec& scalars) +void +ScattDataLegendre::combine(std::vector& those_scatts, + double_1dvec& scalars) { // Find the max order in the data set and make sure we can combine the sets int max_order = 0; @@ -383,8 +396,10 @@ void ScattDataLegendre::combine(std::vector& those_scatts, init(in_gmin, in_gmax, sparse_mult, sparse_scatter); } +//============================================================================== -double_3dvec ScattDataLegendre::get_matrix(int max_order) +double_3dvec +ScattDataLegendre::get_matrix(int max_order) { // Get the sizes and initialize the data to 0 int groups = energy.size(); @@ -408,8 +423,9 @@ double_3dvec ScattDataLegendre::get_matrix(int max_order) // ScattDataHistogram methods //============================================================================== -void ScattDataHistogram::init(int_1dvec& in_gmin, int_1dvec& in_gmax, - double_2dvec& in_mult, double_3dvec& coeffs) +void +ScattDataHistogram::init(int_1dvec& in_gmin, int_1dvec& in_gmax, + double_2dvec& in_mult, double_3dvec& coeffs) { int groups = coeffs.size(); int order = coeffs[0][0].size(); @@ -485,8 +501,10 @@ void ScattDataHistogram::init(int_1dvec& in_gmin, int_1dvec& in_gmax, } } +//============================================================================== -double ScattDataHistogram::calc_f(int gin, int gout, double mu) +double +ScattDataHistogram::calc_f(int gin, int gout, double mu) { double f; if ((gout < gmin[gin]) || (gout > gmax[gin])) { @@ -507,8 +525,10 @@ double ScattDataHistogram::calc_f(int gin, int gout, double mu) return f; } +//============================================================================== -void ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt) +void +ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt) { // Sample the outgoing energy using the base-class method int i_gout; @@ -540,8 +560,10 @@ void ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt) wgt *= mult[gin][i_gout]; } +//============================================================================== -double_3dvec ScattDataHistogram::get_matrix(int max_order) +double_3dvec +ScattDataHistogram::get_matrix(int max_order) { // Get the sizes and initialize the data to 0 int groups = energy.size(); @@ -562,9 +584,11 @@ double_3dvec ScattDataHistogram::get_matrix(int max_order) return matrix; } +//============================================================================== -void ScattDataHistogram::combine(std::vector& those_scatts, - double_1dvec& scalars) +void +ScattDataHistogram::combine(std::vector& those_scatts, + double_1dvec& scalars) { // Find the max order in the data set and make sure we can combine the sets int max_order; @@ -704,8 +728,9 @@ void ScattDataHistogram::combine(std::vector& those_scatts, // ScattDataTabular methods //============================================================================== -void ScattDataTabular::init(int_1dvec& in_gmin, int_1dvec& in_gmax, - double_2dvec& in_mult, double_3dvec& coeffs) +void +ScattDataTabular::init(int_1dvec& in_gmin, int_1dvec& in_gmax, + double_2dvec& in_mult, double_3dvec& coeffs) { int groups = coeffs.size(); int order = coeffs[0][0].size(); @@ -790,8 +815,10 @@ void ScattDataTabular::init(int_1dvec& in_gmin, int_1dvec& in_gmax, } } +//============================================================================== -double ScattDataTabular::calc_f(int gin, int gout, double mu) +double +ScattDataTabular::calc_f(int gin, int gout, double mu) { double f; if ((gout < gmin[gin]) || (gout > gmax[gin])) { @@ -813,8 +840,10 @@ double ScattDataTabular::calc_f(int gin, int gout, double mu) return f; } +//============================================================================== -void ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt) +void +ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt) { // Sample the outgoing energy using the base-class method int i_gout; @@ -859,8 +888,10 @@ void ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt) wgt *= mult[gin][i_gout]; } +//============================================================================== -double_3dvec ScattDataTabular::get_matrix(int max_order) +double_3dvec +ScattDataTabular::get_matrix(int max_order) { // Get the sizes and initialize the data to 0 int groups = energy.size(); @@ -881,8 +912,11 @@ double_3dvec ScattDataTabular::get_matrix(int max_order) return matrix; } -void ScattDataTabular::combine(std::vector& those_scatts, - double_1dvec& scalars) +//============================================================================== + +void +ScattDataTabular::combine(std::vector& those_scatts, + double_1dvec& scalars) { // Find the max order in the data set and make sure we can combine the sets int max_order; @@ -1018,9 +1052,13 @@ void ScattDataTabular::combine(std::vector& those_scatts, init(in_gmin, in_gmax, sparse_mult, sparse_scatter); } +//============================================================================== +// module-level methods +//============================================================================== -void convert_legendre_to_tabular(ScattDataLegendre& leg, - ScattDataTabular& tab, int n_mu) +void +convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab, + int n_mu) { tab.generic_init(n_mu, leg.gmin, leg.gmax, leg.energy, leg.mult); tab.scattxs = leg.scattxs; diff --git a/src/scattdata.h b/src/scattdata.h index 456610655f..e5f115e615 100644 --- a/src/scattdata.h +++ b/src/scattdata.h @@ -47,6 +47,10 @@ class ScattData { virtual double_3dvec get_matrix(int max_order) = 0; }; +//============================================================================== +// ScattDataLegendre represents the angular distributions as Legendre kernels +//============================================================================== + class ScattDataLegendre: public ScattData { protected: // Maximal value for rejection sampling from a rectangle @@ -64,6 +68,11 @@ class ScattDataLegendre: public ScattData { double_3dvec get_matrix(int max_order); }; +//============================================================================== +// ScattDataHistogram represents the angular distributions as a histogram, as it +// would be if it came from a "mu" tally in OpenMC +//============================================================================== + class ScattDataHistogram: public ScattData { protected: double_1dvec mu; @@ -79,6 +88,11 @@ class ScattDataHistogram: public ScattData { double_3dvec get_matrix(int max_order); }; +//============================================================================== +// ScattDataTabular represents the angular distributions as a table of mu and +// f(mu) +//============================================================================== + class ScattDataTabular: public ScattData { protected: double_1dvec mu; @@ -96,8 +110,13 @@ class ScattDataTabular: public ScattData { double_3dvec get_matrix(int max_order); }; -void convert_legendre_to_tabular(ScattDataLegendre& leg, - ScattDataTabular& tab, int n_mu); +//============================================================================== +// Function to convert Legendre functions to tabular +//============================================================================== + +void +convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab, + int n_mu); } // namespace openmc #endif // SCATTDATA_H \ No newline at end of file diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 8fe5e2da1f..806bcc65e8 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -72,9 +72,6 @@ module simulation_header logical :: trace - ! Number of distribcell maps - integer :: n_maps - !$omp threadprivate(trace, thread_id, current_work) contains diff --git a/src/string_utils.h b/src/string_utils.h new file mode 100644 index 0000000000..e2cb773ab1 --- /dev/null +++ b/src/string_utils.h @@ -0,0 +1,35 @@ +#ifndef STRING_UTILS_H +#define STRING_UTILS_H + +#include +#include + + +namespace openmc { + +std::vector +split(const std::string &in) +{ + std::vector out; + + for (int i = 0; i < in.size(); ) { + // Increment i until we find a non-whitespace character. + if (std::isspace(in[i])) { + i++; + + } else { + // Find the next whitespace character at j. + int j = i + 1; + while (j < in.size() && std::isspace(in[j]) == 0) {j++;} + + // Push-back everything between i and j. + out.push_back(in.substr(i, j-i)); + i = j + 1; // j is whitespace so leapfrog to j+1 + } + } + + return out; +} + +} // namespace openmc +#endif // STRING_UTILS_H diff --git a/src/summary.F90 b/src/summary.F90 index cae81a88bc..062b4efc24 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -15,7 +15,6 @@ module summary use surface_header use string, only: to_str use tally_header, only: TallyObject - use tally_filter_distribcell, only: find_offset implicit none private @@ -162,8 +161,7 @@ contains subroutine write_geometry(file_id) integer(HID_T), intent(in) :: file_id - integer :: i, j, k, m - integer, allocatable :: lattice_universes(:,:,:) + integer :: i, j, cell_fill integer, allocatable :: cell_materials(:) integer, allocatable :: cell_ids(:) real(8), allocatable :: cell_temperatures(:) @@ -171,8 +169,7 @@ contains integer(HID_T) :: cells_group, cell_group integer(HID_T) :: surfaces_group integer(HID_T) :: universes_group, univ_group - integer(HID_T) :: lattices_group, lattice_group - character(:), allocatable :: region_spec + integer(HID_T) :: lattices_group type(Cell), pointer :: c class(Lattice), pointer :: lat @@ -181,7 +178,7 @@ contains call write_attribute(geom_group, "n_cells", n_cells) call write_attribute(geom_group, "n_surfaces", n_surfaces) call write_attribute(geom_group, "n_universes", n_universes) - call write_attribute(geom_group, "n_lattices", n_lattices) + call write_attribute(geom_group, "n_lattices", size(lattices)) ! ========================================================================== ! WRITE INFORMATION ON CELLS @@ -192,16 +189,12 @@ contains ! Write information on each cell CELL_LOOP: do i = 1, n_cells c => cells(i) - cell_group = create_group(cells_group, "cell " // trim(to_str(c%id))) + cell_group = create_group(cells_group, "cell " // trim(to_str(c%id()))) - ! Write name for this cell - call write_dataset(cell_group, "name", c%name) - - ! Write universe for this cell - call write_dataset(cell_group, "universe", universes(c%universe)%id) + call c % to_hdf5(cell_group) ! Write information on what fills this cell - select case (c%type) + select case (c%type()) case (FILL_MATERIAL) call write_dataset(cell_group, "fill_type", "material") @@ -233,12 +226,7 @@ contains case (FILL_UNIVERSE) call write_dataset(cell_group, "fill_type", "universe") - call write_dataset(cell_group, "fill", universes(c%fill)%id) - if (allocated(c%offset)) then - if (size(c%offset) > 0) then - call write_dataset(cell_group, "offset", c%offset) - end if - end if + call write_dataset(cell_group, "fill", universes(c%fill()+1)%id) if (allocated(c%translation)) then call write_dataset(cell_group, "translation", c%translation) @@ -249,30 +237,12 @@ contains case (FILL_LATTICE) call write_dataset(cell_group, "fill_type", "lattice") - call write_dataset(cell_group, "lattice", lattices(c%fill)%obj%id) + ! Do not access the 'lattices' array with 'c % fill() + 1' directly; it + ! causes a segfault in GCC 7.3.0 + cell_fill = c % fill() + 1 + call write_dataset(cell_group, "lattice", lattices(cell_fill)%obj%id()) end select - ! Write list of bounding surfaces - region_spec = "" - do j = 1, size(c%region) - k = c%region(j) - select case(k) - case (OP_LEFT_PAREN) - region_spec = trim(region_spec) // " (" - case (OP_RIGHT_PAREN) - region_spec = trim(region_spec) // " )" - case (OP_COMPLEMENT) - region_spec = trim(region_spec) // " ~" - case (OP_INTERSECTION) - case (OP_UNION) - region_spec = trim(region_spec) // " |" - case default - region_spec = trim(region_spec) // " " // to_str(& - sign(surfaces(abs(k))%id(), k)) - end select - end do - call write_dataset(cell_group, "region", adjustl(region_spec)) - call close_group(cell_group) end do CELL_LOOP @@ -307,7 +277,7 @@ contains if (size(u % cells) > 0) then allocate(cell_ids(size(u % cells))) do j = 1, size(u % cells) - cell_ids(j) = cells(u % cells(j)) % id + cell_ids(j) = cells(u % cells(j)) % id() end do call write_dataset(univ_group, "cells", cell_ids) deallocate(cell_ids) @@ -322,87 +292,12 @@ contains ! ========================================================================== ! WRITE INFORMATION ON LATTICES - ! Create lattices group (nothing directly written here) then close lattices_group = create_group(geom_group, "lattices") - ! Write information on each lattice - LATTICE_LOOP: do i = 1, n_lattices + do i = 1, size(lattices) lat => lattices(i)%obj - lattice_group = create_group(lattices_group, "lattice " // trim(to_str(lat%id))) - - ! Write name, pitch, and outer universe - call write_dataset(lattice_group, "name", lat%name) - call write_dataset(lattice_group, "pitch", lat%pitch) - if (lat % outer > 0) then - call write_dataset(lattice_group, "outer", universes(lat % outer) % id) - else - call write_dataset(lattice_group, "outer", lat % outer) - end if - - select type (lat) - type is (RectLattice) - ! Write lattice type. - call write_dataset(lattice_group, "type", "rectangular") - - ! Write lattice dimensions, lower left corner, and pitch - if (lat % is_3d) then - call write_dataset(lattice_group, "dimension", lat % n_cells) - call write_dataset(lattice_group, "lower_left", lat % lower_left) - else - call write_dataset(lattice_group, "dimension", lat % n_cells(1:2)) - call write_dataset(lattice_group, "lower_left", lat % lower_left) - end if - - ! Write lattice universes. - allocate(lattice_universes(lat%n_cells(1), lat%n_cells(2), & - &lat%n_cells(3))) - do j = 1, lat%n_cells(1) - do k = 0, lat%n_cells(2) - 1 - do m = 1, lat%n_cells(3) - lattice_universes(j, k+1, m) = & - universes(lat%universes(j, lat%n_cells(2) - k, m))%id - end do - end do - end do - - type is (HexLattice) - ! Write lattice type. - call write_dataset(lattice_group, "type", "hexagonal") - - ! Write number of lattice cells. - call write_dataset(lattice_group, "n_rings", lat%n_rings) - call write_dataset(lattice_group, "n_axial", lat%n_axial) - - ! Write lattice center - call write_dataset(lattice_group, "center", lat%center) - - ! Write lattice universes. - allocate(lattice_universes(2*lat%n_rings - 1, 2*lat%n_rings - 1, & - &lat%n_axial)) - do m = 1, lat%n_axial - do k = 1, 2*lat%n_rings - 1 - do j = 1, 2*lat%n_rings - 1 - if (j + k < lat%n_rings + 1) then - ! This array position is never used; put a -1 to indicate this - lattice_universes(j,k,m) = -1 - cycle - else if (j + k > 3*lat%n_rings - 1) then - ! This array position is never used; put a -1 to indicate this - lattice_universes(j,k,m) = -1 - cycle - end if - lattice_universes(j,k,m) = universes(lat%universes(j,k,m))%id - end do - end do - end do - end select - - ! Write lattice universes - call write_dataset(lattice_group, "universes", lattice_universes) - deallocate(lattice_universes) - - call close_group(lattice_group) - end do LATTICE_LOOP + call lat % to_hdf5(lattices_group) + end do call close_group(lattices_group) call close_group(geom_group) @@ -504,7 +399,7 @@ contains call write_dataset(material_group, "nuclides", nuc_names) ! Deallocate temporary array deallocate(nuc_names) - ! Write atom densities + ! Write nuclide atom densities call write_dataset(material_group, "nuclide_densities", nuc_densities) deallocate(nuc_densities) end if diff --git a/src/surface.cpp b/src/surface.cpp index 2c3279604a..56e57a854f 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -11,12 +11,25 @@ namespace openmc { +//============================================================================== +// Module constant definitions +//============================================================================== + +extern "C" const int BC_TRANSMIT {0}; +extern "C" const int BC_VACUUM {1}; +extern "C" const int BC_REFLECT {2}; +extern "C" const int BC_PERIODIC {3}; + //============================================================================== // Global variables //============================================================================== int32_t n_surfaces; +Surface **surfaces_c; + +std::map surface_map; + //============================================================================== // Helper functions for reading the "coeffs" node of an XML surface element //============================================================================== @@ -305,7 +318,7 @@ void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "x-plane", false); std::array coeffs {{x0}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } bool SurfaceXPlane::periodic_translate(PeriodicSurface *other, double xyz[3], @@ -370,7 +383,7 @@ void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "y-plane", false); std::array coeffs {{y0}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } bool SurfaceYPlane::periodic_translate(PeriodicSurface *other, double xyz[3], @@ -436,7 +449,7 @@ void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "z-plane", false); std::array coeffs {{z0}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } bool SurfaceZPlane::periodic_translate(PeriodicSurface *other, double xyz[3], @@ -497,7 +510,7 @@ void SurfacePlane::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "plane", false); std::array coeffs {{A, B, C, D}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } bool SurfacePlane::periodic_translate(PeriodicSurface *other, double xyz[3], @@ -629,7 +642,7 @@ void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "x-cylinder", false); std::array coeffs {{y0, z0, r}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } //============================================================================== @@ -663,7 +676,7 @@ void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "y-cylinder", false); std::array coeffs {{x0, z0, r}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } //============================================================================== @@ -697,7 +710,7 @@ void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "z-cylinder", false); std::array coeffs {{x0, y0, r}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } //============================================================================== @@ -768,7 +781,7 @@ void SurfaceSphere::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "sphere", false); std::array coeffs {{x0, y0, z0, r}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } //============================================================================== @@ -885,7 +898,7 @@ void SurfaceXCone::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "x-cone", false); std::array coeffs {{x0, y0, z0, r_sq}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } //============================================================================== @@ -919,7 +932,7 @@ void SurfaceYCone::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "y-cone", false); std::array coeffs {{x0, y0, z0, r_sq}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } //============================================================================== @@ -953,7 +966,7 @@ void SurfaceZCone::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "z-cone", false); std::array coeffs {{x0, y0, z0, r_sq}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } //============================================================================== @@ -1047,7 +1060,7 @@ void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "quadric", false); std::array coeffs {{A, B, C, D, E, F, G, H, J, K}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } //============================================================================== @@ -1116,12 +1129,12 @@ read_surfaces(pugi::xml_node *node) } } - // Fill the surface dictionary. + // Fill the surface map. for (int i_surf = 0; i_surf < n_surfaces; i_surf++) { int id = surfaces_c[i_surf]->id; - auto in_dict = surface_dict.find(id); - if (in_dict == surface_dict.end()) { - surface_dict[id] = i_surf; + auto in_map = surface_map.find(id); + if (in_map == surface_map.end()) { + surface_map[id] = i_surf; } else { std::stringstream err_msg; err_msg << "Two or more surfaces use the same unique ID: " << id; @@ -1211,7 +1224,7 @@ read_surfaces(pugi::xml_node *node) } } else { // Convert the surface id to an index. - surf->i_periodic = surface_dict[surf->i_periodic]; + surf->i_periodic = surface_map[surf->i_periodic]; } } else { // This is a SurfacePlane. We won't try to find it's partner if the @@ -1223,7 +1236,7 @@ read_surfaces(pugi::xml_node *node) fatal_error(err_msg); } else { // Convert the surface id to an index. - surf->i_periodic = surface_dict[surf->i_periodic]; + surf->i_periodic = surface_map[surf->i_periodic]; } } @@ -1242,43 +1255,36 @@ read_surfaces(pugi::xml_node *node) // Fortran compatibility functions //============================================================================== -extern "C" Surface* surface_pointer(int surf_ind) {return surfaces_c[surf_ind];} +extern "C" { + Surface* surface_pointer(int surf_ind) {return surfaces_c[surf_ind];} -extern "C" int surface_id(Surface *surf) {return surf->id;} + int surface_id(Surface *surf) {return surf->id;} -extern "C" int surface_bc(Surface *surf) {return surf->bc;} + int surface_bc(Surface *surf) {return surf->bc;} -extern "C" bool surface_sense(Surface *surf, double xyz[3], double uvw[3]) -{return surf->sense(xyz, uvw);} + void surface_reflect(Surface *surf, double xyz[3], double uvw[3]) + {surf->reflect(xyz, uvw);} -extern "C" void surface_reflect(Surface *surf, double xyz[3], double uvw[3]) -{surf->reflect(xyz, uvw);} + void surface_normal(Surface *surf, double xyz[3], double uvw[3]) + {return surf->normal(xyz, uvw);} -extern "C" double -surface_distance(Surface *surf, double xyz[3], double uvw[3], bool coincident) -{return surf->distance(xyz, uvw, coincident);} + void surface_to_hdf5(Surface *surf, hid_t group) {surf->to_hdf5(group);} -extern "C" void surface_normal(Surface *surf, double xyz[3], double uvw[3]) -{return surf->normal(xyz, uvw);} + int surface_i_periodic(PeriodicSurface *surf) {return surf->i_periodic;} -extern "C" void surface_to_hdf5(Surface *surf, hid_t group) -{surf->to_hdf5(group);} + bool + surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, double xyz[3], + double uvw[3]) + {return surf->periodic_translate(other, xyz, uvw);} -extern "C" int surface_i_periodic(PeriodicSurface *surf) -{return surf->i_periodic;} - -extern "C" bool -surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, double xyz[3], - double uvw[3]) -{return surf->periodic_translate(other, xyz, uvw);} - -extern "C" void free_memory_surfaces_c() -{ - for (int i = 0; i < n_surfaces; i++) {delete surfaces_c[i];} - delete surfaces_c; - surfaces_c = nullptr; - n_surfaces = 0; - surface_dict.clear(); + void free_memory_surfaces_c() + { + for (int i = 0; i < n_surfaces; i++) {delete surfaces_c[i];} + delete surfaces_c; + surfaces_c = nullptr; + n_surfaces = 0; + surface_map.clear(); + } } } // namespace openmc diff --git a/src/surface.h b/src/surface.h index 15ac0d5bdd..1d8cc6ac7b 100644 --- a/src/surface.h +++ b/src/surface.h @@ -8,25 +8,19 @@ #include "hdf5.h" #include "pugixml/pugixml.hpp" +#include "constants.h" + namespace openmc { //============================================================================== -// Module constants +// Module constant declarations (defined in .cpp) //============================================================================== -extern "C" const int BC_TRANSMIT {0}; -extern "C" const int BC_VACUUM {1}; -extern "C" const int BC_REFLECT {2}; -extern "C" const int BC_PERIODIC {3}; - -//============================================================================== -// Constants that should eventually be moved out of this file -//============================================================================== - -extern "C" double FP_COINCIDENT; -constexpr double INFTY {std::numeric_limits::max()}; -constexpr int C_NONE {-1}; +extern "C" const int BC_TRANSMIT; +extern "C" const int BC_VACUUM; +extern "C" const int BC_REFLECT; +extern "C" const int BC_PERIODIC; //============================================================================== // Global variables @@ -35,9 +29,9 @@ constexpr int C_NONE {-1}; extern "C" int32_t n_surfaces; class Surface; -Surface **surfaces_c; +extern Surface **surfaces_c; -std::map surface_dict; +extern std::map surface_map; //============================================================================== //! Coordinates for an axis-aligned cube that bounds a geometric object. @@ -64,7 +58,7 @@ public: //int neighbor_pos[], //!< List of cells on positive side // neighbor_neg[]; //!< List of cells on negative side int bc; //!< Boundary condition - std::string name{""}; //!< User-defined name + std::string name; //!< User-defined name explicit Surface(pugi::xml_node surf_node); @@ -107,6 +101,7 @@ public: //! Write all information needed to reconstruct the surface to an HDF5 group. //! @param group_id An HDF5 group id. + //TODO: this probably needs to include i_periodic for PeriodicSurface void to_hdf5(hid_t group_id) const; protected: @@ -383,19 +378,21 @@ public: // Fortran compatibility functions //============================================================================== -extern "C" Surface* surface_pointer(int surf_ind); -extern "C" int surface_id(Surface *surf); -extern "C" int surface_bc(Surface *surf); -extern "C" bool surface_sense(Surface *surf, double xyz[3], double uvw[3]); -extern "C" void surface_reflect(Surface *surf, double xyz[3], double uvw[3]); -extern "C" double surface_distance(Surface *surf, double xyz[3], double uvw[3], - bool coincident); -extern "C" void surface_normal(Surface *surf, double xyz[3], double uvw[3]); -extern "C" void surface_to_hdf5(Surface *surf, hid_t group); -extern "C" int surface_i_periodic(PeriodicSurface *surf); -extern "C" bool surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, - double xyz[3], double uvw[3]); -extern "C" void free_memory_surfaces_c(); +extern "C" { + Surface* surface_pointer(int surf_ind); + int surface_id(Surface *surf); + int surface_bc(Surface *surf); + bool surface_sense(Surface *surf, double xyz[3], double uvw[3]); + void surface_reflect(Surface *surf, double xyz[3], double uvw[3]); + double surface_distance(Surface *surf, double xyz[3], double uvw[3], + bool coincident); + void surface_normal(Surface *surf, double xyz[3], double uvw[3]); + void surface_to_hdf5(Surface *surf, hid_t group); + int surface_i_periodic(PeriodicSurface *surf); + bool surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, + double xyz[3], double uvw[3]); + void free_memory_surfaces_c(); +} } // namespace openmc #endif // SURFACE_H diff --git a/src/surface_header.F90 b/src/surface_header.F90 index ab6babb5cb..916f9c11c9 100644 --- a/src/surface_header.F90 +++ b/src/surface_header.F90 @@ -30,16 +30,6 @@ module surface_header integer(C_INT) :: bc end function surface_bc_c - pure function surface_sense_c(surf_ptr, xyz, uvw) & - bind(C, name='surface_sense') result(sense) - use ISO_C_BINDING - implicit none - type(C_PTR), intent(in), value :: surf_ptr - real(C_DOUBLE), intent(in) :: xyz(3) - real(C_DOUBLE), intent(in) :: uvw(3) - logical(C_BOOL) :: sense - end function surface_sense_c - pure subroutine surface_reflect_c(surf_ptr, xyz, uvw) & bind(C, name='surface_reflect') use ISO_C_BINDING @@ -49,17 +39,6 @@ module surface_header real(C_DOUBLE), intent(inout) :: uvw(3); end subroutine surface_reflect_c - pure function surface_distance_c(surf_ptr, xyz, uvw, coincident) & - bind(C, name='surface_distance') result(d) - use ISO_C_BINDING - implicit none - type(C_PTR), intent(in), value :: surf_ptr - real(C_DOUBLE), intent(in) :: xyz(3); - real(C_DOUBLE), intent(in) :: uvw(3); - logical(C_BOOL), intent(in), value :: coincident; - real(C_DOUBLE) :: d; - end function surface_distance_c - pure subroutine surface_normal_c(surf_ptr, xyz, uvw) & bind(C, name='surface_normal') use ISO_C_BINDING @@ -115,9 +94,7 @@ module surface_header procedure :: id => surface_id procedure :: bc => surface_bc - procedure :: sense => surface_sense procedure :: reflect => surface_reflect - procedure :: distance => surface_distance procedure :: normal => surface_normal procedure :: to_hdf5 => surface_to_hdf5 procedure :: i_periodic => surface_i_periodic @@ -146,14 +123,6 @@ contains bc = surface_bc_c(this % ptr) end function surface_bc - pure function surface_sense(this, xyz, uvw) result(sense) - class(Surface), intent(in) :: this - real(C_DOUBLE), intent(in) :: xyz(3) - real(C_DOUBLE), intent(in) :: uvw(3) - logical(C_BOOL) :: sense - sense = surface_sense_c(this % ptr, xyz, uvw) - end function surface_sense - pure subroutine surface_reflect(this, xyz, uvw) class(Surface), intent(in) :: this real(C_DOUBLE), intent(in) :: xyz(3); @@ -161,15 +130,6 @@ contains call surface_reflect_c(this % ptr, xyz, uvw) end subroutine surface_reflect - pure function surface_distance(this, xyz, uvw, coincident) result(d) - class(Surface), intent(in) :: this - real(C_DOUBLE), intent(in) :: xyz(3); - real(C_DOUBLE), intent(in) :: uvw(3); - logical(C_BOOL), intent(in) :: coincident; - real(C_DOUBLE) :: d; - d = surface_distance_c(this % ptr, xyz, uvw, coincident) - end function surface_distance - pure subroutine surface_normal(this, xyz, uvw) class(Surface), intent(in) :: this real(C_DOUBLE), intent(in) :: xyz(3); diff --git a/src/tallies/tally_filter_cell.F90 b/src/tallies/tally_filter_cell.F90 index 00ab8fa15b..be37467073 100644 --- a/src/tallies/tally_filter_cell.F90 +++ b/src/tallies/tally_filter_cell.F90 @@ -14,13 +14,14 @@ module tally_filter_cell implicit none private + public :: openmc_cell_filter_get_bins !=============================================================================== ! CELLFILTER specifies which geometric cells tally events reside in. !=============================================================================== type, public, extends(TallyFilter) :: CellFilter - integer, allocatable :: cells(:) + integer(C_INT32_T), allocatable :: cells(:) type(DictIntInt) :: map contains procedure :: from_xml @@ -79,7 +80,7 @@ contains allocate(cell_ids(size(this % cells))) do i = 1, size(this % cells) - cell_ids(i) = cells(this % cells(i)) % id + cell_ids(i) = cells(this % cells(i)) % id() end do call write_dataset(filter_group, "bins", cell_ids) end subroutine to_statepoint_cell @@ -113,7 +114,31 @@ contains integer, intent(in) :: bin character(MAX_LINE_LEN) :: label - label = "Cell " // to_str(cells(this % cells(bin)) % id) + label = "Cell " // to_str(cells(this % cells(bin)) % id()) end function text_label_cell +!=============================================================================== +! C API FUNCTIONS +!=============================================================================== + + function openmc_cell_filter_get_bins(index, cells, n) result(err) bind(C) + ! Return the cells associated with a cell filter + integer(C_INT32_T), value :: index + type(C_PTR), intent(out) :: cells + integer(C_INT32_T), intent(out) :: n + integer(C_INT) :: err + + err = verify_filter(index) + if (err == 0) then + select type (f => filters(index) % obj) + type is (CellFilter) + cells = C_LOC(f % cells) + n = size(f % cells) + class default + err = E_INVALID_TYPE + call set_errmsg("Tried to get cells from a non-cell filter.") + end select + end if + end function openmc_cell_filter_get_bins + end module tally_filter_cell diff --git a/src/tallies/tally_filter_cellborn.F90 b/src/tallies/tally_filter_cellborn.F90 index 0d7461ef82..a373c2d4e1 100644 --- a/src/tallies/tally_filter_cellborn.F90 +++ b/src/tallies/tally_filter_cellborn.F90 @@ -74,7 +74,7 @@ contains call write_dataset(filter_group, "n_bins", this % n_bins) allocate(cell_ids(size(this % cells))) do i = 1, size(this % cells) - cell_ids(i) = cells(this % cells(i)) % id + cell_ids(i) = cells(this % cells(i)) % id() end do call write_dataset(filter_group, "bins", cell_ids) end subroutine to_statepoint_cellborn @@ -108,7 +108,7 @@ contains integer, intent(in) :: bin character(MAX_LINE_LEN) :: label - label = "Birth Cell " // to_str(cells(this % cells(bin)) % id) + label = "Birth Cell " // to_str(cells(this % cells(bin)) % id()) end function text_label_cellborn end module tally_filter_cellborn diff --git a/src/tallies/tally_filter_cellfrom.F90 b/src/tallies/tally_filter_cellfrom.F90 index 2c54021226..83a3179d26 100644 --- a/src/tallies/tally_filter_cellfrom.F90 +++ b/src/tallies/tally_filter_cellfrom.F90 @@ -64,7 +64,7 @@ contains allocate(cell_ids(size(this % cells))) do i = 1, size(this % cells) - cell_ids(i) = cells(this % cells(i)) % id + cell_ids(i) = cells(this % cells(i)) % id() end do call write_dataset(filter_group, "bins", cell_ids) end subroutine to_statepoint_cell_from @@ -74,7 +74,7 @@ contains integer, intent(in) :: bin character(MAX_LINE_LEN) :: label - label = "Cell from " // to_str(cells(this % cells(bin)) % id) + label = "Cell from " // to_str(cells(this % cells(bin)) % id()) end function text_label_cell_from end module tally_filter_cellfrom diff --git a/src/tallies/tally_filter_distribcell.F90 b/src/tallies/tally_filter_distribcell.F90 index 6c42161ca7..c9564083da 100644 --- a/src/tallies/tally_filter_distribcell.F90 +++ b/src/tallies/tally_filter_distribcell.F90 @@ -14,7 +14,6 @@ module tally_filter_distribcell implicit none private - public :: find_offset !=============================================================================== ! DISTRIBCELLFILTER specifies which distributed geometric cells tally events @@ -58,20 +57,20 @@ contains distribcell_index = cells(this % cell) % distribcell_index offset = 0 do i = 1, p % n_coord - if (cells(p % coord(i) % cell) % type == FILL_UNIVERSE) then - offset = offset + cells(p % coord(i) % cell) % & - offset(distribcell_index) - elseif (cells(p % coord(i) % cell) % type == FILL_LATTICE) then + if (cells(p % coord(i) % cell) % type() == FILL_UNIVERSE) then + offset = offset + cells(p % coord(i) % cell) & + % offset(distribcell_index-1) + elseif (cells(p % coord(i) % cell) % type() == FILL_LATTICE) then if (lattices(p % coord(i + 1) % lattice) % obj & % are_valid_indices([& p % coord(i + 1) % lattice_x, & p % coord(i + 1) % lattice_y, & p % coord(i + 1) % lattice_z])) then - offset = offset + lattices(p % coord(i + 1) % lattice) % obj % & - offset(distribcell_index, & - p % coord(i + 1) % lattice_x, & + offset = offset + lattices(p % coord(i + 1) % lattice) % obj & + % offset(distribcell_index - 1, & + [p % coord(i + 1) % lattice_x, & p % coord(i + 1) % lattice_y, & - p % coord(i + 1) % lattice_z) + p % coord(i + 1) % lattice_z]) end if end if if (this % cell == p % coord(i) % cell) then @@ -88,7 +87,7 @@ contains call write_dataset(filter_group, "type", "distribcell") call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", cells(this % cell) % id) + call write_dataset(filter_group, "bins", cells(this % cell) % id()) end subroutine to_statepoint_distribcell subroutine initialize_distribcell(this) @@ -102,7 +101,7 @@ contains val = cell_dict % get(id) if (val /= EMPTY) then this % cell = val - this % n_bins = cells(this % cell) % instances + this % n_bins = cells(this % cell) % n_instances() else call fatal_error("Could not find cell " // trim(to_str(id)) & &// " specified on tally filter.") @@ -114,13 +113,8 @@ contains integer, intent(in) :: bin character(MAX_LINE_LEN) :: label - integer :: offset - type(Universe), pointer :: univ - - univ => universes(root_universe) - offset = 0 label = '' - call find_offset(this % cell, univ, bin-1, offset, label) + call find_offset(this % cell, bin-1, label) label = "Distributed Cell " // label end function text_label_distribcell @@ -130,279 +124,46 @@ contains ! the target cell with the given offset !=============================================================================== - recursive subroutine find_offset(i_cell, univ, target_offset, offset, path) - + subroutine find_offset(i_cell, target_offset, path) integer, intent(in) :: i_cell ! The target cell index - type(Universe), intent(in) :: univ ! Universe to begin search - integer, intent(in) :: target_offset ! Target offset - integer, intent(inout) :: offset ! Current offset - character(*), intent(inout) :: path ! Path to offset + integer, intent(in) :: target_offset ! Target offset + character(*), intent(inout) :: path ! Path to offset integer :: map ! Index in maps vector - integer :: i, j ! Index over cells - integer :: k, l, m ! Indices in lattice - integer :: old_k, old_l, old_m ! Previous indices in lattice - integer :: n_x, n_y, n_z ! Lattice cell array dimensions - integer :: n ! Number of cells to search - integer :: cell_index ! Index in cells array - integer :: lat_offset ! Offset from lattice - integer :: temp_offset ! Looped sum of offsets - integer :: i_univ ! index in universes array - logical :: this_cell = .false. ! Advance in this cell? - logical :: later_cell = .false. ! Fill cells after this one? - type(Cell), pointer :: c ! Pointer to current cell - type(Universe), pointer :: next_univ ! Next universe to loop through - class(Lattice), pointer :: lat ! Pointer to current lattice + integer :: i ! Index over cells + + integer(C_INT) :: path_len + character(kind=C_CHAR), allocatable, target :: path_c(:) + + interface + function distribcell_path_len(target_cell, map, target_offset, root_univ)& + bind(C) result(len) + import C_INT32_T, C_INT + integer(C_INT32_T), intent(in), value :: target_cell, map, & + target_offset, root_univ + integer(C_INT) :: len + end function distribcell_path_len + + subroutine distribcell_path(target_cell, map, target_offset, root_univ, & + path) bind(C) + import C_INT32_T, C_CHAR + integer(C_INT32_T), intent(in), value :: target_cell, map, & + target_offset, root_univ + character(kind=C_CHAR), intent(out) :: path(*) + end subroutine distribcell_path + end interface ! Get the distribcell index for this cell map = cells(i_cell) % distribcell_index - n = size(univ % cells) - - ! Write to the geometry stack - i_univ = universe_dict % get(univ % id) - if (i_univ == root_universe) then - path = trim(path) // "u" // to_str(univ%id) - else - path = trim(path) // "->u" // to_str(univ%id) - end if - - ! Look through all cells in this universe - do i = 1, n - ! If the cell matches the goal and the offset matches final, write to the - ! geometry stack - if (univ % cells(i) == i_cell .and. offset == target_offset) then - c => cells(univ % cells(i)) - path = trim(path) // "->c" // to_str(c % id) - return - end if - end do - - ! Find the fill cell or lattice cell that we need to enter - do i = 1, n - - later_cell = .false. - - c => cells(univ % cells(i)) - - this_cell = .false. - - ! If we got here, we still think the target is in this universe - ! or further down, but it's not this exact cell. - ! Compare offset to next cell to see if we should enter this cell - if (i /= n) then - - do j = i+1, n - - c => cells(univ % cells(j)) - - ! Skip normal cells which do not have offsets - if (c % type == FILL_MATERIAL) cycle - - ! Break loop once we've found the next cell with an offset - exit - end do - - ! Ensure we didn't just end the loop by iteration - if (c % type /= FILL_MATERIAL) then - - ! There are more cells in this universe that it could be in - later_cell = .true. - - ! Two cases, lattice or fill cell - if (c % type == FILL_UNIVERSE) then - temp_offset = c % offset(map) - - ! Get the offset of the first lattice location - else - lat => lattices(c % fill) % obj - temp_offset = lat % offset(map, 1, 1, 1) - end if - - ! If the final offset is in the range of offset - temp_offset+offset - ! then the goal is in this cell - if (target_offset < temp_offset + offset) then - this_cell = .true. - end if - end if - end if - - if (n == 1 .and. c % type /= FILL_MATERIAL) then - this_cell = .true. - end if - - if (.not. later_cell) then - this_cell = .true. - end if - - ! Get pointer to THIS cell because target must be in this cell - if (this_cell) then - - cell_index = univ % cells(i) - c => cells(cell_index) - - path = trim(path) // "->c" // to_str(c%id) - - ! ==================================================================== - ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL - if (c % type == FILL_UNIVERSE) then - - ! Enter this cell to update the current offset - offset = c % offset(map) + offset - - next_univ => universes(c % fill) - call find_offset(i_cell, next_univ, target_offset, offset, path) - return - - ! ==================================================================== - ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL - elseif (c % type == FILL_LATTICE) then - - ! Set current lattice - lat => lattices(c % fill) % obj - - select type (lat) - - ! ================================================================== - ! RECTANGULAR LATTICES - type is (RectLattice) - - ! Write to the geometry stack - path = trim(path) // "->l" // to_str(lat%id) - - n_x = lat % n_cells(1) - n_y = lat % n_cells(2) - n_z = lat % n_cells(3) - old_m = 1 - old_l = 1 - old_k = 1 - - ! Loop over lattice coordinates - do k = 1, n_x - do l = 1, n_y - do m = 1, n_z - - if (target_offset >= lat % offset(map, k, l, m) + offset) then - if (k == n_x .and. l == n_y .and. m == n_z) then - ! This is last lattice cell, so target must be here - lat_offset = lat % offset(map, k, l, m) - offset = offset + lat_offset - next_univ => universes(lat % universes(k, l, m)) - if (lat % is_3d) then - path = trim(path) // "(" // trim(to_str(k-1)) // & - "," // trim(to_str(l-1)) // "," // & - trim(to_str(m-1)) // ")" - else - path = trim(path) // "(" // trim(to_str(k-1)) // & - "," // trim(to_str(l-1)) // ")" - end if - call find_offset(i_cell, next_univ, target_offset, offset, path) - return - else - old_m = m - old_l = l - old_k = k - cycle - end if - else - ! Target is at this lattice position - lat_offset = lat % offset(map, old_k, old_l, old_m) - offset = offset + lat_offset - next_univ => universes(lat % universes(old_k, old_l, old_m)) - if (lat % is_3d) then - path = trim(path) // "(" // trim(to_str(old_k-1)) // & - "," // trim(to_str(old_l-1)) // "," // & - trim(to_str(old_m-1)) // ")" - else - path = trim(path) // "(" // trim(to_str(old_k-1)) // & - "," // trim(to_str(old_l-1)) // ")" - end if - call find_offset(i_cell, next_univ, target_offset, offset, path) - return - end if - - end do - end do - end do - - ! ================================================================== - ! HEXAGONAL LATTICES - type is (HexLattice) - - ! Write to the geometry stack - path = trim(path) // "->l" // to_str(lat%id) - - n_z = lat % n_axial - n_y = 2 * lat % n_rings - 1 - n_x = 2 * lat % n_rings - 1 - old_m = 1 - old_l = 1 - old_k = 1 - - ! Loop over lattice coordinates - do m = 1, n_z - do l = 1, n_y - do k = 1, n_x - - ! This array position is never used - if (k + l < lat % n_rings + 1) then - cycle - ! This array position is never used - else if (k + l > 3*lat % n_rings - 1) then - cycle - end if - - if (target_offset >= lat % offset(map, k, l, m) + offset) then - if (k == lat % n_rings .and. l == n_y .and. m == n_z) then - ! This is last lattice cell, so target must be here - lat_offset = lat % offset(map, k, l, m) - offset = offset + lat_offset - next_univ => universes(lat % universes(k, l, m)) - if (lat % is_3d) then - path = trim(path) // "(" // & - trim(to_str(k - lat % n_rings)) // "," // & - trim(to_str(l - lat % n_rings)) // "," // & - trim(to_str(m - 1)) // ")" - else - path = trim(path) // "(" // & - trim(to_str(k - lat % n_rings)) // "," // & - trim(to_str(l - lat % n_rings)) // ")" - end if - call find_offset(i_cell, next_univ, target_offset, offset, path) - return - else - old_m = m - old_l = l - old_k = k - cycle - end if - else - ! Target is at this lattice position - lat_offset = lat % offset(map, old_k, old_l, old_m) - offset = offset + lat_offset - next_univ => universes(lat % universes(old_k, old_l, old_m)) - if (lat % is_3d) then - path = trim(path) // "(" // & - trim(to_str(old_k - lat % n_rings)) // "," // & - trim(to_str(old_l - lat % n_rings)) // "," // & - trim(to_str(old_m - 1)) // ")" - else - path = trim(path) // "(" // & - trim(to_str(old_k - lat % n_rings)) // "," // & - trim(to_str(old_l - lat % n_rings)) // ")" - end if - call find_offset(i_cell, next_univ, target_offset, offset, path) - return - end if - - end do - end do - end do - - end select - - end if - end if + path_len = distribcell_path_len(i_cell-1, map-1, target_offset, & + root_universe-1) + allocate(path_c(path_len)) + call distribcell_path(i_cell-1, map-1, target_offset, root_universe-1, & + path_c) + do i = 1, min(path_len, MAX_LINE_LEN) + if (path_c(i) == C_NULL_CHAR) exit + path(i:i) = path_c(i) end do end subroutine find_offset diff --git a/src/tracking.F90 b/src/tracking.F90 index 80fb55a800..c926739848 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -188,7 +188,7 @@ contains p % coord(p % n_coord) % cell = NONE if (any(lattice_translation /= 0)) then ! Particle crosses lattice boundary - p % surface = NONE + p % surface = ERROR_INT call cross_lattice(p, lattice_translation) p % event = EVENT_LATTICE else @@ -219,7 +219,7 @@ contains call score_surface_tally(p, active_meshsurf_tallies) ! Clear surface component - p % surface = NONE + p % surface = ERROR_INT if (run_CE) then call collision(p) @@ -489,7 +489,7 @@ contains ! COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS ! Remove lower coordinate levels and assignment of surface - p % surface = NONE + p % surface = ERROR_INT p % n_coord = 1 call find_cell(p, found) diff --git a/src/volume_calc.F90 b/src/volume_calc.F90 index 8264797a7a..f793d4689e 100644 --- a/src/volume_calc.F90 +++ b/src/volume_calc.F90 @@ -205,7 +205,8 @@ contains elseif (this % domain_type == FILTER_CELL) THEN do level = 1, p % n_coord do i_domain = 1, size(this % domain_id) - if (cells(p % coord(level) % cell) % id == this % domain_id(i_domain)) then + if (cells(p % coord(level) % cell) % id() & + == this % domain_id(i_domain)) then i_material = p % material call check_hit(i_domain, i_material, indices, hits, n_mat) end if diff --git a/src/xml_interface.cpp b/src/xml_interface.cpp index aea1c2c17b..b4fcc524c6 100644 --- a/src/xml_interface.cpp +++ b/src/xml_interface.cpp @@ -2,9 +2,7 @@ #include // for std::transform #include -#include -#include "pugixml/pugixml.hpp" #include "error.h" diff --git a/src/xml_interface.h b/src/xml_interface.h index 52aa6003f4..b6f9e3246d 100644 --- a/src/xml_interface.h +++ b/src/xml_interface.h @@ -2,6 +2,7 @@ #define XML_INTERFACE_H #include +#include #include "pugixml/pugixml.hpp" diff --git a/src/xsdata.cpp b/src/xsdata.cpp index 655bf9a891..ddfe4c0fdf 100644 --- a/src/xsdata.cpp +++ b/src/xsdata.cpp @@ -66,10 +66,12 @@ XsData::XsData(int energy_groups, int num_delayed_groups, bool fissionable, } } +//============================================================================== -void XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format, - int final_scatter_format, int order_data, int max_order, - int legendre_to_tabular_points, bool is_isotropic) +void +XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format, + int final_scatter_format, int order_data, int max_order, + int legendre_to_tabular_points, bool is_isotropic) { // Reconstruct the dimension information so it doesn't need to be passed int n_pol = total.size(); @@ -127,8 +129,10 @@ void XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format, } } +//============================================================================== -void XsData::_fissionable_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, +void +XsData::_fissionable_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, int energy_groups, int delayed_groups, bool is_isotropic) { @@ -525,8 +529,10 @@ void XsData::_fissionable_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, } } +//============================================================================== -void XsData::_scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, +void +XsData::_scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, int energy_groups, int scatter_format, int final_scatter_format, int order_data, int max_order, int legendre_to_tabular_points) { @@ -656,14 +662,16 @@ void XsData::_scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, for (int p = 0; p < n_pol; p++) { for (int a = 0; a < n_azi; a++) { scatter[p][a]->init(gmin[p][a], gmax[p][a], temp_mult[p][a], - input_scatt[p][a]); + input_scatt[p][a]); } } } } +//============================================================================== -void XsData::combine(std::vector those_xs, double_1dvec& scalars) +void +XsData::combine(std::vector those_xs, double_1dvec& scalars) { // Combine the non-scattering data for (int i = 0; i < those_xs.size(); i++) { @@ -751,8 +759,10 @@ void XsData::combine(std::vector those_xs, double_1dvec& scalars) } } +//============================================================================== -bool XsData::equiv(const XsData& that) +bool +XsData::equiv(const XsData& that) { bool match = false; // check n_pol (total.size()), n_azi (total[0].size()), and diff --git a/tests/dummy_operator.py b/tests/dummy_operator.py index 05fe97a950..c66070ad40 100644 --- a/tests/dummy_operator.py +++ b/tests/dummy_operator.py @@ -17,8 +17,8 @@ class DummyOperator(TransportOperator): y_2(1.5) ~ 3.1726475740397628 """ - def __init__(self): - pass + def __init__(self, previous_results=None): + self.prev_res = previous_results def __call__(self, vec, power, print_out=False): """Evaluates F(y) diff --git a/tests/regression_tests/asymmetric_lattice/results_true.dat b/tests/regression_tests/asymmetric_lattice/results_true.dat index ef7608f2fc..f60cfa64f5 100644 --- a/tests/regression_tests/asymmetric_lattice/results_true.dat +++ b/tests/regression_tests/asymmetric_lattice/results_true.dat @@ -1 +1 @@ -3b76b468b9c0e7df6508189b75748a1b7b4f2b37486396749e1a15e536526336f70b60bb207095c39ecbd46822e8c8705ea81184a3c8546e7da09bb905d74637 \ No newline at end of file +bd3fd10177bc0c7b8542f15228decf8608de013872d201e6ae885cf9b5b6d7516ddb32044ee2f7fd5b207e7d37d5ce3f03e347a2e850953bba4fc26b150c89d2 \ No newline at end of file diff --git a/tests/regression_tests/distribmat/inputs_true.dat b/tests/regression_tests/distribmat/inputs_true.dat index 39e611e16a..a1c1e83116 100644 --- a/tests/regression_tests/distribmat/inputs_true.dat +++ b/tests/regression_tests/distribmat/inputs_true.dat @@ -1,7 +1,7 @@ - + diff --git a/tests/regression_tests/distribmat/results_true.dat b/tests/regression_tests/distribmat/results_true.dat index cd6aa95d4f..e0143f0fff 100644 --- a/tests/regression_tests/distribmat/results_true.dat +++ b/tests/regression_tests/distribmat/results_true.dat @@ -3,7 +3,7 @@ k-combined: Cell ID = 11 Name = - Fill = [2, 3, None, 2] + Fill = [2, None, 3, 2] Region = -9 Rotation = None Translation = None diff --git a/tests/regression_tests/distribmat/test.py b/tests/regression_tests/distribmat/test.py index de1b778772..cb598b9d3f 100644 --- a/tests/regression_tests/distribmat/test.py +++ b/tests/regression_tests/distribmat/test.py @@ -34,7 +34,7 @@ class DistribmatTestHarness(PyAPITestHarness): r0 = openmc.ZCylinder(R=0.3) c11 = openmc.Cell(cell_id=11, region=-r0) - c11.fill = [dense_fuel, light_fuel, None, dense_fuel] + c11.fill = [dense_fuel, None, light_fuel, dense_fuel] c12 = openmc.Cell(cell_id=12, region=+r0, fill=moderator) fuel_univ = openmc.Universe(universe_id=11, cells=[c11, c12]) diff --git a/tests/regression_tests/filter_distribcell/case-1/results_true.dat b/tests/regression_tests/filter_distribcell/case-1/results_true.dat index a824d7fce3..2d9835b1ab 100644 --- a/tests/regression_tests/filter_distribcell/case-1/results_true.dat +++ b/tests/regression_tests/filter_distribcell/case-1/results_true.dat @@ -3,10 +3,10 @@ k-combined: tally 1: 1.548980E-02 2.399339E-04 -1.278781E-02 -1.635280E-04 1.426319E-02 2.034385E-04 +1.278781E-02 +1.635280E-04 1.018927E-02 1.038213E-04 tally 2: diff --git a/tests/regression_tests/filter_distribcell/case-3/results_true.dat b/tests/regression_tests/filter_distribcell/case-3/results_true.dat index 8eb7c2b609..acd74781db 100644 --- a/tests/regression_tests/filter_distribcell/case-3/results_true.dat +++ b/tests/regression_tests/filter_distribcell/case-3/results_true.dat @@ -1 +1 @@ -386147796cf64c908002a81bc47af6763a14811ba6c457ca790653ceb8ef1c6b0376b94783f9746baf7389f18321b788c72758ec56ffd05fa3146139e6e53308 \ No newline at end of file +a1f6ccf0bef1075ffc12e7fa058a44c83360801ba6e2c76f00e1bd7e9543fe1832777b122321dcb1cfb5fc1165ce84f711a65ca83903dd57cd09652e003daf62 \ No newline at end of file diff --git a/tests/regression_tests/multipole/inputs_true.dat b/tests/regression_tests/multipole/inputs_true.dat index d1ef3e00aa..3bbcc20247 100644 --- a/tests/regression_tests/multipole/inputs_true.dat +++ b/tests/regression_tests/multipole/inputs_true.dat @@ -1,7 +1,7 @@ - + diff --git a/tests/regression_tests/multipole/results_true.dat b/tests/regression_tests/multipole/results_true.dat index 0ab4191d8b..d734e67bd0 100644 --- a/tests/regression_tests/multipole/results_true.dat +++ b/tests/regression_tests/multipole/results_true.dat @@ -37,5 +37,5 @@ Cell Fill = Material 2 Region = -1 Rotation = None - Temperature = [ 500. 0. 700. 800.] + Temperature = [ 500. 700. 0. 800.] Translation = None diff --git a/tests/regression_tests/multipole/test.py b/tests/regression_tests/multipole/test.py index 5c79d4d6b6..ef797a1bc5 100644 --- a/tests/regression_tests/multipole/test.py +++ b/tests/regression_tests/multipole/test.py @@ -29,7 +29,7 @@ def make_model(): r0 = openmc.ZCylinder(R=0.3) c11 = openmc.Cell(cell_id=11, fill=dense_fuel, region=-r0) - c11.temperature = [500, 0, 700, 800] + c11.temperature = [500, 700, 0, 800] c12 = openmc.Cell(cell_id=12, fill=moderator, region=+r0) fuel_univ = openmc.Universe(universe_id=11, cells=(c11, c12)) diff --git a/tests/regression_tests/tally_aggregation/results_true.dat b/tests/regression_tests/tally_aggregation/results_true.dat index ee763b3bac..0f0a3d60a9 100644 --- a/tests/regression_tests/tally_aggregation/results_true.dat +++ b/tests/regression_tests/tally_aggregation/results_true.dat @@ -1 +1 @@ -eb2002dd2f3016154e7014501629227eb2e5b2a9655b5c0cb3b9d4d681f918fae4197bf4756fe9865c8d34f392b981b287a15e9b2fab5a46b2a5b8a33a8ae770 \ No newline at end of file +8294c7481b3433a4beb25541ae72c1ea915def0c0d0f393f7585371877187f4a2a25e70bb602232e2bc1e877e78db0e9384591e5c71575659f95abb10fae0af3 \ No newline at end of file diff --git a/tests/regression_tests/test_deplete_full.py b/tests/regression_tests/test_deplete_full.py index 19e2c7664b..58e6ddac40 100644 --- a/tests/regression_tests/test_deplete_full.py +++ b/tests/regression_tests/test_deplete_full.py @@ -22,6 +22,7 @@ def test_full(run_in_tmpdir): This test runs a complete OpenMC simulation and tests the outputs. It will take a while. + """ n_rings = 2 diff --git a/tests/regression_tests/test_reference.h5 b/tests/regression_tests/test_reference.h5 index 5323a9a904..10dc37b511 100644 Binary files a/tests/regression_tests/test_reference.h5 and b/tests/regression_tests/test_reference.h5 differ diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index a21e858fda..af013cbb5b 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -26,6 +26,15 @@ def pincell_model(): mat_tally.scores = ['total', 'elastic', '(n,gamma)'] pincell.tallies.append(mat_tally) + # Add an expansion tally + zernike_tally = openmc.Tally() + filter3 = openmc.ZernikeFilter(5, r=.63) + cells = pincell.geometry.root_universe.cells + filter4 = openmc.CellFilter(list(cells.values())) + zernike_tally.filters = [filter3, filter4] + zernike_tally.scores = ['fission'] + pincell.tallies.append(zernike_tally) + # Write XML files in tmpdir with cdtemp(): pincell.export_to_xml() @@ -132,7 +141,7 @@ def test_settings(capi_init): def test_tally_mapping(capi_init): tallies = openmc.capi.tallies assert isinstance(tallies, Mapping) - assert len(tallies) == 1 + assert len(tallies) == 2 for tally_id, tally in tallies.items(): assert isinstance(tally, openmc.capi.Tally) assert tally_id == tally.id @@ -168,6 +177,14 @@ def test_tally(capi_init): t.active = True assert t.active + t2 = openmc.capi.tallies[2] + t2.id = 2 + assert len(t2.filters) == 2 + assert isinstance(t2.filters[0], openmc.capi.ZernikeFilter) + assert isinstance(t2.filters[1], openmc.capi.CellFilter) + assert len(t2.filters[1].bins) == 3 + assert t2.filters[0].order == 5 + def test_new_tally(capi_init): with pytest.raises(exc.AllocationError): @@ -176,7 +193,7 @@ def test_new_tally(capi_init): new_tally.scores = ['flux'] new_tally_with_id = openmc.capi.Tally(10) new_tally_with_id.scores = ['flux'] - assert len(openmc.capi.tallies) == 3 + assert len(openmc.capi.tallies) == 4 def test_tally_results(capi_run): @@ -187,6 +204,10 @@ def test_tally_results(capi_run): assert np.all(t.std_dev[nonzero] >= 0) assert np.all(t.ci_width()[nonzero] >= 1.95*t.std_dev[nonzero]) + t2 = openmc.capi.tallies[2] + n = 5 + assert t2.mean.size == (n + 1) * (n + 2) // 2 * 3 # Number of Zernike coeffs * 3 cells + def test_global_tallies(capi_run): assert openmc.capi.num_realizations() == 5 diff --git a/tests/unit_tests/test_deplete_integrator.py b/tests/unit_tests/test_deplete_integrator.py index a1768625b2..8cade45b14 100644 --- a/tests/unit_tests/test_deplete_integrator.py +++ b/tests/unit_tests/test_deplete_integrator.py @@ -25,6 +25,9 @@ def test_results_save(run_in_tmpdir): # Mock geometry op = MagicMock() + # Avoid DummyOperator thinking it's doing a restart calculation + op.prev_res = None + vol_dict = {} full_burn_list = [] @@ -72,8 +75,8 @@ def test_results_save(run_in_tmpdir): op_result1 = [OperatorResult(k, rates) for k, rates in zip(eigvl1, rate1)] op_result2 = [OperatorResult(k, rates) for k, rates in zip(eigvl2, rate2)] - Results.save(op, x1, op_result1, t1, 0) - Results.save(op, x2, op_result2, t2, 1) + Results.save(op, x1, op_result1, t1, 0, 0) + Results.save(op, x2, op_result2, t2, 0, 1) # Load the files res = ResultsList("depletion_results.h5") diff --git a/tests/unit_tests/test_deplete_predictor.py b/tests/unit_tests/test_deplete_predictor.py index 50803e5085..b39cc7dcac 100644 --- a/tests/unit_tests/test_deplete_predictor.py +++ b/tests/unit_tests/test_deplete_predictor.py @@ -10,7 +10,7 @@ from tests import dummy_operator def test_predictor(run_in_tmpdir): - """Integral regression test of integrator algorithm using predictor/corrector""" + """Integral regression test of integrator algorithm using predictor""" op = dummy_operator.DummyOperator() op.output_dir = "test_integrator_regression" diff --git a/tests/unit_tests/test_deplete_restart.py b/tests/unit_tests/test_deplete_restart.py new file mode 100644 index 0000000000..745f2cce5b --- /dev/null +++ b/tests/unit_tests/test_deplete_restart.py @@ -0,0 +1,168 @@ +"""Regression tests for openmc.deplete restart capability. + +These tests run in two steps, a first run then a restart run, a simple test +problem described in dummy_geometry.py. +""" + +from pytest import approx +import openmc.deplete + +from tests import dummy_operator + + +def test_restart_predictor(run_in_tmpdir): + """Integral regression test of integrator algorithm using predictor.""" + + op = dummy_operator.DummyOperator() + output_dir = "test_restart_predictor" + op.output_dir = output_dir + + # Perform simulation using the predictor algorithm + dt = [0.75] + power = 1.0 + openmc.deplete.predictor(op, dt, power, print_out=False) + + # Load the files + prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + # Re-create depletion operator and load previous results + op = dummy_operator.DummyOperator(prev_res) + op.output_dir = output_dir + + # Perform restarts simulation using the predictor algorithm + openmc.deplete.predictor(op, dt, power, print_out=False) + + # Load the files + res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + _, y1 = res.get_atoms("1", "1") + _, y2 = res.get_atoms("1", "2") + + # Mathematica solution + s1 = [2.46847546272295, 0.986431226850467] + s2 = [4.11525874568034, -0.0581692232513460] + + assert y1[1] == approx(s1[0]) + assert y2[1] == approx(s1[1]) + + assert y1[2] == approx(s2[0]) + assert y2[2] == approx(s2[1]) + + +def test_restart_cecm(run_in_tmpdir): + """Integral regression test of integrator algorithm using CE/CM.""" + + op = dummy_operator.DummyOperator() + output_dir = "test_restart_cecm" + op.output_dir = output_dir + + # Perform simulation using the MCNPX/MCNP6 algorithm + dt = [0.75] + power = 1.0 + openmc.deplete.cecm(op, dt, power, print_out=False) + + # Load the files + prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + # Re-create depletion operator and load previous results + op = dummy_operator.DummyOperator(prev_res) + op.output_dir = output_dir + + # Perform restarts simulation using the MCNPX/MCNP6 algorithm + openmc.deplete.cecm(op, dt, power, print_out=False) + + # Load the files + res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + _, y1 = res.get_atoms("1", "1") + _, y2 = res.get_atoms("1", "2") + + # Mathematica solution + s1 = [1.86872629872102, 1.395525772416039] + s2 = [2.18097439443550, 2.69429754646747] + + assert y1[1] == approx(s1[0]) + assert y2[1] == approx(s1[1]) + + assert y1[3] == approx(s2[0]) + assert y2[3] == approx(s2[1]) + + +def test_restart_predictor_cecm(run_in_tmpdir): + """Integral regression test of integrator algorithm using predictor + for the first run then CE/CM for the restart run.""" + + op = dummy_operator.DummyOperator() + output_dir = "test_restart_predictor_cecm" + op.output_dir = output_dir + + # Perform simulation using the predictor algorithm + dt = [0.75] + power = 1.0 + openmc.deplete.predictor(op, dt, power, print_out=False) + + # Load the files + prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + # Re-create depletion operator and load previous results + op = dummy_operator.DummyOperator(prev_res) + op.output_dir = output_dir + + # Perform restarts simulation using the MCNPX/MCNP6 algorithm + openmc.deplete.cecm(op, dt, power, print_out=False) + + # Load the files + res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + _, y1 = res.get_atoms("1", "1") + _, y2 = res.get_atoms("1", "2") + + # Test solution + s1 = [2.46847546272295, 0.986431226850467] + s2 = [3.09106948392, 0.607102912398] + + assert y1[1] == approx(s1[0]) + assert y2[1] == approx(s1[1]) + + assert y1[2] == approx(s2[0]) + assert y2[2] == approx(s2[1]) + + +def test_restart_cecm_predictor(run_in_tmpdir): + """Integral regression test of integrator algorithm using CE/CM for the + first run then predictor for the restart run.""" + + op = dummy_operator.DummyOperator() + output_dir = "test_restart_cecm_predictor" + op.output_dir = output_dir + + # Perform simulation using the MCNPX/MCNP6 algorithm + dt = [0.75] + power = 1.0 + openmc.deplete.cecm(op, dt, power, print_out=False) + + # Load the files + prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + # Re-create depletion operator and load previous results + op = dummy_operator.DummyOperator(prev_res) + op.output_dir = output_dir + + # Perform restarts simulation using the predictor algorithm + openmc.deplete.predictor(op, dt, power, print_out=False) + + # Load the files + res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + _, y1 = res.get_atoms("1", "1") + _, y2 = res.get_atoms("1", "2") + + # Test solution + s1 = [1.86872629872102, 1.395525772416039] + s2 = [3.32776806576, 2.391425905] + + assert y1[1] == approx(s1[0]) + assert y2[1] == approx(s1[1]) + + assert y1[2] == approx(s2[0]) + assert y2[2] == approx(s2[1]) diff --git a/tests/unit_tests/test_deplete_resultslist.py b/tests/unit_tests/test_deplete_resultslist.py index aad8cd9f68..86896cc151 100644 --- a/tests/unit_tests/test_deplete_resultslist.py +++ b/tests/unit_tests/test_deplete_resultslist.py @@ -22,8 +22,8 @@ def test_get_atoms(res): n_ref = [6.6747328233649218e+08, 3.5421791038348462e+14, 3.6208592242443462e+14, 3.3799758969347038e+14] - np.testing.assert_array_equal(t, t_ref) - np.testing.assert_array_equal(n, n_ref) + np.testing.assert_allclose(t, t_ref) + np.testing.assert_allclose(n, n_ref) def test_get_reaction_rate(res): """Tests evaluating reaction rate.""" @@ -35,8 +35,8 @@ def test_get_reaction_rate(res): xs_ref = np.array([4.0594392323131994e-05, 3.9249546927524987e-05, 3.8394587728581798e-05, 4.1521845978371697e-05]) - np.testing.assert_array_equal(t, t_ref) - np.testing.assert_array_equal(r, n_ref * xs_ref) + np.testing.assert_allclose(t, t_ref) + np.testing.assert_allclose(r, n_ref * xs_ref) def test_get_eigenvalue(res): @@ -47,5 +47,5 @@ def test_get_eigenvalue(res): k_ref = [1.181281798790367, 1.1798750921988739, 1.1965943696058159, 1.2207119847790813] - np.testing.assert_array_equal(t, t_ref) - np.testing.assert_array_equal(k, k_ref) + np.testing.assert_allclose(t, t_ref) + np.testing.assert_allclose(k, k_ref) diff --git a/tests/unit_tests/test_mesh_from_lattice.py b/tests/unit_tests/test_mesh_from_lattice.py new file mode 100644 index 0000000000..4d79fc579d --- /dev/null +++ b/tests/unit_tests/test_mesh_from_lattice.py @@ -0,0 +1,116 @@ +import numpy as np +import openmc +import pytest + + +@pytest.fixture(scope='module') +def pincell1(uo2, water): + cyl = openmc.ZCylinder(R=0.35) + fuel = openmc.Cell(fill=uo2, region=-cyl) + moderator = openmc.Cell(fill=water, region=+cyl) + + univ = openmc.Universe(cells=[fuel, moderator]) + univ.fuel = fuel + univ.moderator = moderator + return univ + + +@pytest.fixture(scope='module') +def pincell2(uo2, water): + cyl = openmc.ZCylinder(R=0.4) + fuel = openmc.Cell(fill=uo2, region=-cyl) + moderator = openmc.Cell(fill=water, region=+cyl) + + univ = openmc.Universe(cells=[fuel, moderator]) + univ.fuel = fuel + univ.moderator = moderator + return univ + + +@pytest.fixture(scope='module') +def zr(): + zr = openmc.Material() + zr.add_element('Zr', 1.0) + zr.set_density('g/cm3', 1.0) + return zr + + +@pytest.fixture(scope='module') +def rlat2(pincell1, pincell2, uo2, water, zr): + """2D Rectangular lattice for testing.""" + all_zr = openmc.Cell(fill=zr) + pitch = 1.2 + n = 3 + u1, u2 = pincell1, pincell2 + lattice = openmc.RectLattice() + lattice.lower_left = (-pitch*n/2, -pitch*n/2) + lattice.pitch = (pitch, pitch) + lattice.outer = openmc.Universe(cells=[all_zr]) + lattice.universes = [ + [u1, u2, u1], + [u2, u1, u2], + [u2, u1, u1] + ] + + return lattice + + +@pytest.fixture(scope='module') +def rlat3(pincell1, pincell2, uo2, water, zr): + """3D Rectangular lattice for testing.""" + + # Create another universe for top layer + hydrogen = openmc.Material() + hydrogen.add_element('H', 1.0) + hydrogen.set_density('g/cm3', 0.09) + h_cell = openmc.Cell(fill=hydrogen) + u3 = openmc.Universe(cells=[h_cell]) + + all_zr = openmc.Cell(fill=zr) + pitch = 1.2 + n = 3 + u1, u2 = pincell1, pincell2 + lattice = openmc.RectLattice() + lattice.lower_left = (-pitch*n/2, -pitch*n/2, -10.0) + lattice.pitch = (pitch, pitch, 10.0) + lattice.outer = openmc.Universe(cells=[all_zr]) + lattice.universes = [ + [[u1, u2, u1], + [u2, u1, u2], + [u2, u1, u1]], + [[u3, u1, u2], + [u1, u3, u2], + [u2, u1, u1]] + ] + + return lattice + + +def test_mesh2d(rlat2): + shape = np.array(rlat2.shape) + width = shape*rlat2.pitch + + mesh1 = openmc.Mesh.from_rect_lattice(rlat2) + assert np.array_equal(mesh1.dimension, (3, 3)) + assert np.array_equal(mesh1.lower_left, rlat2.lower_left) + assert np.array_equal(mesh1.upper_right, rlat2.lower_left + width) + + mesh2 = openmc.Mesh.from_rect_lattice(rlat2, division=3) + assert np.array_equal(mesh2.dimension, (9, 9)) + assert np.array_equal(mesh2.lower_left, rlat2.lower_left) + assert np.array_equal(mesh2.upper_right, rlat2.lower_left + width) + + +def test_mesh3d(rlat3): + shape = np.array(rlat3.shape) + width = shape*rlat3.pitch + + mesh1 = openmc.Mesh.from_rect_lattice(rlat3) + assert np.array_equal(mesh1.dimension, (3, 3, 2)) + assert np.array_equal(mesh1.lower_left, rlat3.lower_left) + assert np.array_equal(mesh1.upper_right, rlat3.lower_left + width) + + mesh2 = openmc.Mesh.from_rect_lattice(rlat3, division=3) + assert np.array_equal(mesh2.dimension, (9, 9, 6)) + assert np.array_equal(mesh2.lower_left, rlat3.lower_left) + assert np.array_equal(mesh2.upper_right, rlat3.lower_left + width) diff --git a/tests/unit_tests/test_transfer_volumes.py b/tests/unit_tests/test_transfer_volumes.py new file mode 100644 index 0000000000..9128ed9694 --- /dev/null +++ b/tests/unit_tests/test_transfer_volumes.py @@ -0,0 +1,45 @@ +"""Regression tests for openmc.deplete.Results.transfer_volumes method. + + +""" + +from pytest import approx +import openmc.deplete + +from tests import dummy_operator + + +def test_transfer_volumes(run_in_tmpdir): + """Unit test of volume transfer in restart calculations.""" + + op = dummy_operator.DummyOperator() + op.output_dir = "test_transfer_volumes" + + # Perform simulation using the predictor algorithm + dt = [0.75] + power = 1.0 + openmc.deplete.predictor(op, dt, power, print_out=False) + + # Load the files + res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + # Create a dictionary of volumes to transfer + res[0].volume['1'] = 1.5 + res[0].volume['2'] = 2.5 + + # Create dummy geometry + mat1 = openmc.Material(material_id=1) + mat1.depletable = True + mat2 = openmc.Material(material_id=2) + + cell = openmc.Cell() + cell.fill = [mat1, mat2] + root = openmc.Universe() + root.add_cell(cell) + geometry = openmc.Geometry(root) + + # Transfer volumes + res[0].transfer_volumes(geometry) + + assert mat1.volume == 1.5 + assert mat2.volume is None