/
manager.py
365 lines (301 loc) · 12.4 KB
/
manager.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
"""
Depletion manager
Control time-steps, material divisions, depletion, etc
"""
import numbers
from collections.abc import Sequence, Callable
from itertools import repeat, starmap
import multiprocessing
import numpy
from hydep import BurnableMaterial, DepletionChain
from hydep.constants import SECONDS_PER_DAY
from hydep.typed import TypedAttr, IterableOf
from hydep.internal import Cram16Solver, Cram48Solver, CompBundle
from hydep.internal.features import FeatureCollection, MICRO_REACTION_XS, FISSION_YIELDS
from hydep.internal.utils import FakeSequence
__all__ = ["Manager"]
class Manager:
"""Primary depletion manager
Responsible for depleting materials and updating compositions
Parameters
----------
chain : hydep.DepletionChain
Chain describing how isotopes decay and transmute
daysteps : iterable of float
Length in time [d] for each coarse depletion step
power : float or iterable of float
Power [W] to deplete the entire problem. If a single value
is given, a constant power will be used. Otherwise, must
have the same number of elements as ``daysteps``, corresponding
to a constant power in each depletion step
substepDivision : int or sequence of int
Number of substeps to divide each entry in ``daysteps``, also
the number of transport solutions, high fidelity and reduced
order, per entry. If an integer, the value will be applied to
all entries.
numPreliminary : int, optional
Number of coarse depletion steps to take before engaging in
coupled behavior. Useful for approaching some equilibrium value
with smaller steps with the high fidelity code
depletionSolver : string or int or callable, optional
Value to use in configuring the depletion solver. Passed to
:meth:`setDepletionSolver`
Attributes
----------
chain : hydep.DepletionChain
Depletion chain
timesteps : numpy.ndarray
Length of time [s] for each coarse depletion step
powers : numpy.ndarray
Power [W] to use for each coarse depletion steps
numPreliminary : int
Number of coarse depletion steps to take before engaging in
coupled behavior. Useful for approaching some equilibrium value
with smaller steps with the high fidelity code
burnable : None or tuple of hydep.BurnableMaterial
Ordering of burnable materials. Must be set prior to depletion
to maintain a consistent mapping from reaction rates and
compositions
needs : hydep.internal.features.FeatureCollection
Read-only property describing the capabilities any solver must
have in order to properly perform the depletion analysis.
Currently requires calculation of isotopic reaction cross
sections in each burnable material, as well as the flux.
substeps : sequence of int
Number of transport solutions, high fidelity and reduced order,
per coarse depletion step. Not writable.
"""
chain = TypedAttr("chain", DepletionChain)
_burnable = IterableOf("burnable", BurnableMaterial, allowNone=True)
def __init__(
self,
chain,
daysteps,
power,
substepDivision,
numPreliminary=0,
depletionSolver=None,
):
self.chain = chain
daysteps = numpy.asarray(daysteps, dtype=float)
if len(daysteps.shape) > 1:
raise TypeError("Day steps must be vector, not array")
self.timesteps = tuple(daysteps * SECONDS_PER_DAY)
self.powers = tuple(self._validatePowers(power))
self._burnable = None
if numPreliminary is None:
self._nprelim = 0
else:
if not isinstance(numPreliminary, numbers.Integral):
raise TypeError(
f"Non-integer preliminary steps not allowed: {type(numPreliminary)}"
)
elif not (0 <= numPreliminary < len(self.timesteps)):
raise ValueError(
"Number of preliminary steps must be between [0, "
f"{len(self.timesteps)}), not {numPreliminary}"
)
self._nprelim = numPreliminary
self._substeps = self._validateSubsteps(substepDivision)
self.setDepletionSolver(depletionSolver)
def _validatePowers(self, power):
if isinstance(power, numbers.Real):
if power <= 0:
raise ValueError(f"Power must be positive, not {power}")
return repeat(power, len(self.timesteps))
elif isinstance(power, Sequence):
if len(power) != len(self.timesteps):
raise ValueError(
f"Number of powers {len(power)} differ from steps "
f"{len(self.timesteps)}"
)
for p in power:
if not isinstance(p, numbers.Real):
raise TypeError(
"Power must be positive real, or vector of positive real. "
f"Found {p}"
)
elif p <= 0:
raise ValueError(
"Power must be positive real, or vector of positive real. "
f"Found {p}"
)
return power
else:
raise TypeError(
"Power must be positive real, or vector of positive real, "
f"not {type(power)}"
)
def _validateSubsteps(self, divisions):
maxallowed = len(self.timesteps) - self._nprelim
if isinstance(divisions, numbers.Integral):
if divisions <= 0:
raise ValueError(f"Divisions must be positive integer, not {divisions}")
return FakeSequence(divisions, maxallowed)
if isinstance(divisions, (numpy.ndarray, Sequence)):
if len(divisions) != maxallowed:
raise ValueError(
"Number of divisions {} not equal to number of time "
"steps {} - number of preliminary steps {}".format(
len(divisions), len(self.timesteps), self._nprelim
)
)
substeps = []
for value in divisions:
if not isinstance(value, numbers.Integral):
raise TypeError(
f"Divisions must be positive integer, not {type(value)}"
)
elif value <= 0:
raise ValueError(f"Divisions must be positive integer, not {value}")
substeps.append(value)
return tuple(substeps)
raise TypeError(
"Substeps must be postive integer, or sequence of positive "
f"integer, not {divisions}"
)
def setDepletionSolver(self, solver):
"""Configure the depletion solver
Solver can either be a string, e.g. ``"cram16"``,
integer, ``16``, or a callable function. Callable functions
should fulfill the following requirements:
1. Be importable / pickle-able in order to be dispatched via
:meth:`multiprocessing.Pool.starmap`
2. Have a call signature ``solver(A, N0, dt)`` where ``A`` is
the :class:`scipy.sparse.csr_matrix` sparse representation
of the depletion matrix with shape ``N x N``, ``N0`` is a
:class:`numpy.ndarray` with the beginning-of-step
compositions, and ``dt`` is the length of the depletion
interval in seconds
For the time being, no introspection is performed to ensure
that the correct signature is used. String values are
case-insensitive, and integers indicate the order of CRAM
to be used.
Parameters
----------
solver : str or int or callable or None
Item indicating what solver should be used. A value
of ``None`` reverts to the default CRAM16.
Raises
------
TypeError
If ``solver`` doesn't match any requirements
"""
if solver is None:
self._depsolver = Cram16Solver.__call__
return
if isinstance(solver, str):
solver = solver.lower()
candidate = {
"cram16": Cram16Solver,
"cram48": Cram48Solver,
16: Cram16Solver,
48: Cram48Solver,
"16": Cram16Solver,
"48": Cram48Solver,
}.get(solver)
if candidate is not None:
self._depsolver = candidate.__call__
return
if isinstance(solver, Callable):
self._depsolver = solver
return solver
raise TypeError(f"Could not decipher {solver} of type {type(solver)}")
@property
def burnable(self):
return self._burnable
@property
def needs(self):
return FeatureCollection({MICRO_REACTION_XS, FISSION_YIELDS})
@property
def numPreliminary(self):
return self._nprelim
@property
def substeps(self):
return self._substeps
def preliminarySteps(self):
"""Iterate over preliminary time steps and powers
Useful for running only a high fidelity solver before
substep depletion with a reduced order solver
Yields
------
float
Time step [s]
float
Power [W] for current step
"""
return zip(self.timesteps[: self._nprelim], self.powers[: self._nprelim])
def activeSteps(self):
"""Iterate over active time steps and powers
These are steps after the :meth:`preliminarySteps`.
Yields
------
float
Time step [s]
float
Power [W] for current step
"""
return zip(self.timesteps[self._nprelim :], self.powers[self._nprelim :])
def beforeMain(self, model, settings=None):
"""Check that all materials have volumes and set indexes
Parameters
----------
model : hydep.Model
Problem to be solved. Must contain at least one
:class:`hydep.BurnableMaterial`
settings : hydep.settings.HydepSettings, optional
Settings for the framework. Use to configure depletion
solver.
"""
if settings is not None:
solver = settings.depletionSolver
if solver is not None:
self.setDepletionSolver(solver)
burnable = tuple(model.root.findBurnableMaterials())
if not burnable:
raise ValueError(f"No burnable materials found in {model}")
for ix, mat in enumerate(burnable):
if mat.volume is None:
raise AttributeError(
f"{mat.__class__} {mat.name} does not have a volume set"
)
mat.index = ix
self._burnable = burnable
def deplete(self, dtSeconds, concentrations, reactionRates, fissionYields):
"""Deplete all burnable materials
Parameters
----------
dtSeconds : float
Length of depletion interval in seconds
concentrations : hydep.internal.CompBundle
Incoming material compositions. The
:attr:`hydep.internal.CompBundle.isotopes` will be
consistent between the incoming and outgoing bundle
reactionRates : iterable of hydep.internal.MicroXsVector
Stand in for microscopic reaction rates in each burnable
material
fissionYields : iterable of hydep.internal.FissionYield
Fission yields in each burnable material
Returns
-------
hydep.internal.CompBundle
New compositions for each burnable material and the isotope
ordering
"""
nr = len(reactionRates)
nm = len(concentrations.densities)
nf = len(fissionYields)
if not nr == nm == nf:
raise ValueError(
"Inconsistent number of reaction rates {} to burnable "
"materials {} and fission yields {}".format(nr, nm, nf)
)
zaiOrder = {iso.zai: ix for ix, iso in enumerate(concentrations.isotopes)}
matrices = starmap(
self.chain.formMatrix,
zip(reactionRates, fissionYields, repeat(zaiOrder, nm)),
)
inputs = zip(matrices, concentrations.densities, repeat(dtSeconds, nm))
with multiprocessing.Pool() as p:
out = p.starmap(self._depsolver, inputs)
return CompBundle(concentrations.isotopes, out)