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

first commit to new branch

Browse source
This commit is contained in:
church89 2023-08-23 18:03:07 +02:00
parent f14fc55e60
commit 38ed171ceb
1 changed files with 435 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

435
openmc/deplete/batchwise.py Normal file
View file

@ -0,0 +1,435 @@
from abc import ABC, abstractmethod
from collections.abc import Iterable
from openmc.search import _SCALAR_BRACKETED_METHODS, search_for_keff
from openmc import Materials, Material, Cell
from openmc.data import atomic_mass, AVOGADRO, ELEMENT_SYMBOL
import numpy as np
import os
import h5py
class Batchwise(ABC):
def __init__(self, operator, model, bracket, bracket_limit,
bracketed_method='brentq', tol=0.01, target=1.0,
print_iterations=True, search_for_keff_output=True):
self.operator = operator
self.burn_mats = operator.burnable_mats
self.local_mats = operator.local_mats
self.model = model
self.geometry = model.geometry
check_iterable_type('bracket', bracket, Real)
check_length('bracket', bracket, 2)
check_less_than('bracket values', bracket[0], bracket[1])
self.bracket = bracket
check_iterable_type('bracket_limit', bracket_limit, Real)
check_length('bracket_limit', bracket_limit, 2)
check_less_than('bracket limit values',
bracket_limit[0], bracket_limit[1])
self.bracket_limit = bracket_limit
self.bracketed_method = bracketed_method
self.tol = tol
self.target = target
self.print_iterations = print_iterations
self.search_for_keff_output = search_for_keff_output
@property
def bracketed_method(self):
return self._bracketed_method
@bracketed_method.setter
def bracketed_method(self, value):
check_value('bracketed_method', value, _SCALAR_BRACKETED_METHODS)
if value != 'brentq':
warn('brentq bracketed method is recomended here')
self._bracketed_method = value
@property
def tol(self):
return self._tol
@tol.setter
def tol(self, value):
check_type("tol", value, Real)
self._tol = value
@property
def target(self):
return self._target
@target.setter
def target(self, value):
check_type("target", value, Real)
self._target = value
@abstractmethod
def _model_builder(self, param):
def _get_cell_id(self, val):
"""Helper method for getting cell id from cell instance or cell name.
Parameters
----------
val : Openmc.Cell or str or int representing Cell
Returns
-------
id : str
Cell id
"""
if isinstance(val, Cell):
check_value('Cell id', val.id, [cell.id for cell in \
self.geometry.get_all_cells().values()])
val = val.id
elif isinstance(val, str):
if val.isnumeric():
check_value('Cell id', val, [str(cell.id) for cell in \
self.geometry.get_all_cells().values()])
val = int(val)
else:
check_value('Cell name', val, [cell.name for cell in \
self.geometry.get_all_cells().values()])
val = [cell.id for cell in \
self.geometry.get_all_cells().values() \
if cell.name == val][0]
elif isinstance(val, int):
check_value('Cell id', val, [cell.id for cell in \
self.geometry.get_all_cells().values()])
else:
ValueError(f'Cell: {val} is not recognized')
return val
def _search_for_keff(self, val):
"""
Perform the criticality search for a given parametric model.
It supports geometrical or material based `search_for_keff`.
Parameters
----------
val : float
Previous result value
Returns
-------
root : float
Estimated value of the variable parameter where keff is the
targeted value
"""
#make sure we don't modify original bracket and tol values
bracket = deepcopy(self.bracket)
#search_for_keff tolerance should vary according to the first guess value
if abs(val) > 1.0:
tol = self.tol / abs(val)
else:
tol = self.tol
# Run until a search_for_keff root is found or ouf ot limits
root = None
while res == None:
search = search_for_keff(self._model_builder,
bracket = [bracket[0] + val, bracket[1] + val],
tol = tol,
bracketed_method = self.bracketed_method,
target = self.target,
print_iterations = self.print_iterations,
run_args = {'output': self.search_for_keff_output})
# if len(search) is 3 search_for_keff was successful
if len(search) == 3:
res,_,_ = search
#Check if root is within bracket limits
if self.bracket_limit[0] < res < self.bracket_limit[1]:
root = res
else:
# Set res with the closest limit and continue
arg_min = abs(np.array(self.bracket_limit) - res).argmin()
warn('WARNING: Search_for_keff returned root out of '\
'bracket limit. Set root to {:.2f} and continue.'
.format(self.bracket_limit[arg_min]))
root = self.bracket_limit[arg_min]
elif len(search) == 2:
guesses, keffs = search
#Check if all guesses are within bracket limits
if all(self.bracket_limit[0] < guess < self.bracket_limit[1] \
for guess in guesses):
#Simple method to iteratively adapt the bracket
print('INFO: Function returned values below or above ' \
'target. Adapt bracket...')
# if the bracket ends up being smaller than the std of the
# keff's closer value to target, no need to continue-
if all(keff <= max(keffs).s for keff in keffs):
arg_min = abs(self.target - np.array(guesses)).argmin()
root = guesses[arg_min]
# Calculate gradient as ratio of delta bracket and delta keffs
grad = abs(np.diff(bracket) / np.diff(keffs))[0].n
# Move the bracket closer to presumed keff root.
# Two cases: both keffs are below or above target
if np.mean(keffs) < self.target:
# direction of moving bracket: +1 is up, -1 is down
if guess[np.argmax(keffs)] > guess[np.argmin(keffs)]:
dir = 1
else:
dir = -1
bracket[np.argmin(keffs)] = bracket[np.argmax(keffs)]
bracket[np.argmax(keffs)] += grad * (self.target - \
max(keffs).n) * dir
else:
if guess[np.argmax(keffs)] > guess[np.argmin(keffs)]:
dir = -1
else:
dir = 1
bracket[np.argmax(keffs)] = bracket[np.argmin(keffs)]
bracket[np.argmin(keffs)] += grad * (min(keffs).n - \
self.target) * dir
else:
# Set res with closest limit and continue
arg_min = abs(np.array(self.bracket_limit) - guesses).argmin()
warn('WARNING: Adaptive iterative bracket went off '\
'bracket limits. Set root to {:.2f} and continue.'
.format(self.bracket_limit[arg_min]))
root = self.bracket_limit[arg_min]
else:
raise ValueError(f'ERROR: Search_for_keff output is not valid')
return root
def _save_res(self, type, step_index, root):
"""
Save results to msr_results.h5 file.
Parameters
----------
type : str
String to characterize geometry and material results
step_index : int
depletion time step index
root : float or dict
Root of the search_for_keff function
"""
filename = 'msr_results.h5'
kwargs = {'mode': "a" if os.path.isfile(filename) else "w"}
if comm.rank == 0:
with h5py.File(filename, **kwargs) as h5:
name = '_'.join([type, str(step_index)])
if name in list(h5.keys()):
last = sorted([int(re.split('_',i)[1]) for i in h5.keys()])[-1]
step_index = last + 1
h5.create_dataset('_'.join([type, str(step_index)]), data=root)
def _update_volumes_after_depletion(self, x):
"""
After a depletion step, both materials volume and density change, due to
decay, transmutation reactions and transfer rates, if set.
At present we lack an implementation to calculate density and volume
changes due to the different molecules speciation. Therefore, the assumption
is to consider the density constant and let the material volume
vary with the change in nuclides concentrations.
The method uses the nuclides concentrations coming from the previous Bateman
solution and calculates a new volume, keeping the mass density of the material
constant. It will then assign the volumes to the AtomNumber instance.
Parameters
----------
x : list of numpy.ndarray
Total atom concentrations
"""
self.operator.number.set_density(x)
for rank in range(comm.size):
number_i = comm.bcast(self.operator.number, root=rank)
for i, mat in enumerate(number_i.materials):
# Total nuclides density
density = 0
for nuc in number_i.nuclides:
# total number of atoms
val = number_i[mat, nuc]
# obtain nuclide density in atoms-g/mol
density += val * atomic_mass(nuc)
# Get mass dens from beginning, intended to be held constant
rho = openmc.lib.materials[int(mat)].get_density('g/cm3')
number_i.volume[i] = density / AVOGADRO / rho
class BatchwiseCell(Batchwise):
def __init__(self, operator, model, cell, attrib_name, axis, bracket,
bracket_limit, bracketed_method='brentq', tol=0.01, target=1.0,
print_iterations=True, search_for_keff_output=True):
super().__init__(operator, model, bracket, bracket_limit,
bracketed_method, tol, target, print_iterations,
search_for_keff_output)
self.cell_id = super()._get_cell_id(cell_id_or_name)
check_value('attrib_name', attrib_name,
('rotation', 'translation'))
self.attrib_name = attrib_name
#index of cell directionnal axis
check_value('axis', axis, (0,1,2))
self.axis = axis
# Initialize vector
self.vector = np.zeros(3)
# materials that fill the attribute cell, if depletables
self.cell_material_ids = [cell.fill.id for cell in \
self.geometry.get_all_cells()[self.cell_id].fill.cells.values() \
if cell.fill.depletable]
def _get_cell_attrib(self):
"""
Get cell attribute coefficient.
Returns
-------
coeff : float
cell coefficient
"""
for cell in openmc.lib.cells.values():
if cell.id == self.cell_id:
if self.attrib_name == 'translation':
return cell.translation[self.axis]
elif self.attrib_name == 'rotation':
return cell.rotation[self.axis]
def _set_cell_attrib(self, val):
"""
Set cell attribute to the cell instance.
Attributes are only applied to cells filled with a universe
Parameters
----------
var : float
Surface coefficient to set
geometry : openmc.model.geometry
OpenMC geometry model
attrib_name : str
Currently only translation is implemented
"""
self.vector[self.axis] = val
for cell in openmc.lib.cells.values():
if cell.id == self.cell_id or cell.name == self.cell_id:
setattr(cell, self.attrib_name, self.vector)
def _update_materials(self, x):
"""
Assign concentration vectors from Bateman solution at previous
timestep to the in-memory model materials, after having recalculated the
material volume.
Parameters
----------
x : list of numpy.ndarray
Total atom concentrations
"""
super()._update_volumes_after_depletion(x)
for rank in range(comm.size):
number_i = comm.bcast(self.operator.number, root=rank)
for mat in number_i.materials:
nuclides = []
densities = []
for nuc in number_i.nuclides:
# get atom density in atoms/b-cm
val = 1.0e-24 * number_i.get_atom_density(mat, nuc)
if nuc in self.operator.nuclides_with_data:
if val > self.atom_density_limit:
nuclides.append(nuc)
densities.append(val)
#set nuclide densities to model in memory (C-API)
openmc.lib.materials[int(mat)].set_densities(nuclides, densities)
def _update_volumes(self):
openmc.lib.calculate_volumes()
res = openmc.VolumeCalculation.from_hdf5('volume_1.h5')
number_i = self.operator.number
for mat_idx, mat_id in enumerate(self.local_mats):
if mat_id in self.cell_material_ids:
number_i.volume[mat_idx] = res.volumes[int(mat_id)].n
def _update_x(self):
number_i = self.operator.number
for mat_idx, mat_id in enumerate(self.local_mats):
if mat_id in self.cell_material_ids:
for nuc_idx, nuc in enumerate(number_i.burnable_nuclides):
x[mat_idx][nuc_idx] = number_i.volume[mat_idx] * \
number_i.get_atom_density(mat_idx, nuc)
return x
def _model_builder(self, param):
"""
Builds the parametric model that is passed to the `msr_search_for_keff`
function by setting the parametric variable to the geoemetrical cell.
Parameters
----------
param : model parametricl variable
for examlple: cell translation coefficient
Returns
-------
_model : openmc.model.Model
OpenMC parametric model
"""
self._set_cell_attrib(param)
#Calulate new volume and update if materials filling the cell are
# depletable
if self.cell_material_ids:
self._update_volumes()
return self.model
def search_for_keff(self, x, step_index):
"""
Perform the criticality search on the parametric geometrical model.
Will set the root of the `search_for_keff` function to the cell
attribute.
Parameters
----------
x : list of numpy.ndarray
Total atoms concentrations
Returns
-------
x : list of numpy.ndarray
Updated total atoms concentrations
"""
# Get cell attribute from previous iteration
val = self._get_cell_attrib()
check_type('Cell coeff', val, Real)
# Update volume and concentration vectors before performing the search_for_keff
self._update_materials(x)
# Calculate new cell attribute
root = super().search_for_keff(val)
# set results value as attribute in the geometry
self._set_cell_attrib(root)
print('UPDATE: old value: {:.2f} cm --> ' \
'new value: {:.2f} cm'.format(val, root))
# x needs to be updated after the search with new volumes if materials cell
# depletable
if self.cell_material_ids:
self._update_x()
#Store results
super()._save_res('geometry', step_index, root)
return x