Line data Source code
1 : /********************************************************************/
2 : /* SOFTWARE COPYRIGHT NOTIFICATION */
3 : /* Cardinal */
4 : /* */
5 : /* (c) 2021 UChicago Argonne, LLC */
6 : /* ALL RIGHTS RESERVED */
7 : /* */
8 : /* Prepared by UChicago Argonne, LLC */
9 : /* Under Contract No. DE-AC02-06CH11357 */
10 : /* With the U. S. Department of Energy */
11 : /* */
12 : /* Prepared by Battelle Energy Alliance, LLC */
13 : /* Under Contract No. DE-AC07-05ID14517 */
14 : /* With the U. S. Department of Energy */
15 : /* */
16 : /* See LICENSE for full restrictions */
17 : /********************************************************************/
18 :
19 : #pragma once
20 :
21 : #define LIBMESH
22 :
23 : #include "CardinalProblem.h"
24 : #include "PostprocessorInterface.h"
25 : #include "CardinalEnums.h"
26 :
27 : #include "mpi.h"
28 : #include "openmc/bank.h"
29 : #include "openmc/capi.h"
30 : #include "openmc/cell.h"
31 : #include "openmc/geometry.h"
32 : #include "openmc/geometry_aux.h"
33 : #include "openmc/hdf5_interface.h"
34 : #include "openmc/material.h"
35 : #include "openmc/mesh.h"
36 : #include "openmc/settings.h"
37 : #include "openmc/simulation.h"
38 : #include "openmc/source.h"
39 : #include "openmc/state_point.h"
40 : #include "openmc/tallies/tally.h"
41 : #include "openmc/tallies/filter_cell_instance.h"
42 : #include "xtensor/xview.hpp"
43 :
44 : // Forward declarations to avoid cyclic dependencies.
45 : class OpenMCNuclideDensities;
46 : class OpenMCDomainFilterEditor;
47 : class OpenMCTallyEditor;
48 :
49 : /**
50 : * Base class for all MOOSE wrappings of OpenMC
51 : */
52 : class OpenMCProblemBase : public CardinalProblem, public PostprocessorInterface
53 : {
54 : public:
55 : OpenMCProblemBase(const InputParameters & params);
56 :
57 : static InputParameters validParams();
58 :
59 : virtual ~OpenMCProblemBase() override;
60 :
61 : /**
62 : * Get the subdomain name for a given ID. If not named, we return the ID
63 : * @param[in] id subdomain ID
64 : * @return name
65 : */
66 : std::string subdomainName(const SubdomainID & id) const;
67 :
68 : /**
69 : * Print a full error message when catching errors from OpenMC
70 : * @param[in] err OpenMC error code
71 : * @param[in] descriptor descriptive message for error
72 : */
73 : void catchOpenMCError(const int & err, const std::string descriptor) const;
74 :
75 : /**
76 : * Whether the score is a reaction rate score
77 : * @return whether the tally from OpenMC has units of 1/src
78 : */
79 : bool isReactionRateScore(const std::string & score) const;
80 :
81 : /**
82 : * Whether the score is a heating-type score
83 : * @return whether the tally from OpenMC has units of eV/src
84 : */
85 : bool isHeatingScore(const std::string & score) const;
86 :
87 : /**
88 : * Add a constant monomial auxiliary variable
89 : * @param[in] name name of the variable
90 : * @param[in] block optional subdomain names on which to restrict the variable
91 : * @return numeric index for the variable in the auxiliary system
92 : */
93 : unsigned int addExternalVariable(const std::string & name, const std::vector<SubdomainName> * block = nullptr);
94 :
95 : /**
96 : * Get the scaling value applied to the [Mesh] to convert to OpenMC's centimeters units
97 : * @return scaling value
98 : */
99 1291326 : const Real & scaling() const { return _scaling; }
100 :
101 : /**
102 : * Whether the problem has user defined scaling or not.
103 : * @return whether the user has set the problem scaling or not
104 : */
105 340000 : bool hasScaling() const { return _specified_scaling; }
106 :
107 : /**
108 : * Convert from a MOOSE-type enum into a valid OpenMC tally score string
109 : * @param[in] score MOOSE-type enum string
110 : * @return OpenMC tally score string
111 : */
112 : std::string enumToTallyScore(const std::string & score) const;
113 :
114 : /**
115 : * Convert into a MOOSE-type enum from a valid OpenMC tally score string
116 : * @param[in] score OpenMC tally score string
117 : * @return MOOSE-type enum string
118 : */
119 : std::string tallyScoreToEnum(const std::string & score) const;
120 :
121 : /**
122 : * Find the geometry type in the OpenMC model
123 : * @param[out] has_csg_universe whether there is at least 1 CSG universe
124 : * @param[out] has_dag_universe whether there is at least 1 DagMC universe
125 : */
126 : virtual void geometryType(bool & has_csg_universe, bool & has_dag_universe) const;
127 :
128 : /// Whether this is the first time OpenMC is running
129 : bool firstSolve() const;
130 :
131 : /**
132 : * Convert from a MooseEnum for a trigger metric to an OpenMC enum
133 : * @param[in] trigger trigger metric
134 : * @return OpenMC enum
135 : */
136 : openmc::TriggerMetric triggerMetric(trigger::TallyTriggerTypeEnum trigger) const;
137 : openmc::TriggerMetric triggerMetric(std::string trigger) const;
138 :
139 : /**
140 : * Convert from a MooseEnum for tally estimator to an OpenMC enum
141 : * @param[in] estimator MOOSE estimator enum
142 : * @return OpenMC enum
143 : */
144 : openmc::TallyEstimator tallyEstimator(tally::TallyEstimatorEnum estimator) const;
145 :
146 : /**
147 : * Convert a tally estimator to a string (for output purposes).
148 : * @param[in] estimator OpenMC tally estimator enum
149 : * @return a string form of the OpenMC tally estimator enum
150 : */
151 : std::string estimatorToString(openmc::TallyEstimator estimator) const;
152 :
153 : /// Run an OpenMC simulation
154 : void externalSolve() override;
155 :
156 : /// Set up an OpenMC simulation.
157 : virtual void initialSetup() override;
158 :
159 : /// Set the 'mesh changed' adaptivity flag.
160 : virtual void syncSolutions(ExternalProblem::Direction direction) override;
161 : virtual bool adaptMesh() override;
162 :
163 : /// Import temperature and density from a properties.h5 file
164 : void importProperties() const;
165 :
166 : /**
167 : * \brief Compute the sum of a tally within each bin
168 : *
169 : * For example, suppose we have a cell tally with 4 bins, one for each of 4
170 : * different cells. This function will return the sum of the tally in each of
171 : * those bins, so the return xtensor will have a length of 4, with each value
172 : * representing the sum for that bin.
173 : *
174 : * @param[in] tally OpenMC tally
175 : * @param[in] score tally score
176 : * @return tally sum within each bin
177 : */
178 : xt::xtensor<double, 1> tallySum(openmc::Tally * tally, const unsigned int & score) const;
179 :
180 : /**
181 : * Compute the sum of a tally across all of its bins
182 : * @param[in] tally OpenMC tallies (multiple if repeated mesh tallies)
183 : * @param[in] score tally score
184 : * @return tally sum
185 : */
186 : double tallySumAcrossBins(std::vector<openmc::Tally *> tally, const unsigned int & score) const;
187 :
188 : /**
189 : * Compute the mean of a tally across all of its bins
190 : * @param[in] tally OpenMC tallies (multiple if repeated mesh tallies)
191 : * @param[in] score tally score
192 : * @return tally mean
193 : */
194 : double tallyMeanAcrossBins(std::vector<openmc::Tally *> tally, const unsigned int & score) const;
195 :
196 : /**
197 : * Type definition for storing the relevant aspects of the OpenMC geometry; the first
198 : * value is the cell index, while the second is the cell instance.
199 : */
200 : typedef std::pair<int32_t, int32_t> cellInfo;
201 :
202 : /**
203 : * Whether a cell is filled with VOID (vacuum)
204 : * @param[in] cell_info cell ID, instance
205 : * @return whether cell is void
206 : */
207 : bool cellIsVoid(const cellInfo & cell_info) const;
208 :
209 : /**
210 : * Whether this cell has zero instances
211 : * @param[in] cell_info cell info
212 : * @return whether this cell has zero instances
213 : */
214 : bool cellHasZeroInstances(const cellInfo & cell_info) const;
215 :
216 : /**
217 : * Get the material name given its index. If the material does not have a name,
218 : * return the ID.
219 : * @param[in] index
220 : * @return material name
221 : */
222 : std::string materialName(const int32_t index) const;
223 :
224 : /**
225 : * Compute relative error
226 : * @param[in] sum sum of scores
227 : * @param[in] sum_sq sum of scores squared
228 : * @param[in] n_realizations number of realizations
229 : */
230 : xt::xtensor<double, 1> relativeError(const xt::xtensor<double, 1> & sum,
231 : const xt::xtensor<double, 1> & sum_sq, const int & n_realizations) const;
232 :
233 : /**
234 : * Get the density conversion factor (multiplicative factor)
235 : * @return density conversion factor from kg/m3 to g/cm3
236 : */
237 : const Real & densityConversionFactor() const { return _density_conversion_factor; }
238 :
239 : /**
240 : * Get the number of particles used in the current Monte Carlo calculation
241 : * @return number of particles
242 : */
243 : int nParticles() const;
244 :
245 : /**
246 : * Total number of particles run (not multiplied by batches)
247 : * @return total number of particles
248 : */
249 24 : int nTotalParticles() const { return _total_n_particles; }
250 :
251 : /**
252 : * Get the cell ID from the cell index
253 : * @param[in] index cell index
254 : * @return cell ID
255 : */
256 : int32_t cellID(const int32_t index) const;
257 :
258 : /**
259 : * Get the material ID from the material index; for VOID cells, this returns -1
260 : * @param[in] index material index
261 : * @return cell material ID
262 : */
263 : int32_t materialID(const int32_t index) const;
264 :
265 : /**
266 : * Print point coordinates with a neater formatting than the default MOOSE printing
267 : * @return formatted point string
268 : */
269 : std::string printPoint(const Point & p) const;
270 :
271 : /**
272 : * Get a descriptive, formatted, string describing a material
273 : * @param[in] index material index
274 : * @return descriptive string
275 : */
276 : std::string printMaterial(const int32_t & index) const;
277 :
278 : /**
279 : * Write the source bank to HDF5 for postprocessing or for use in subsequent solves
280 : * @param[in] filename file name
281 : */
282 : void writeSourceBank(const std::string & filename);
283 :
284 : /**
285 : * Get the total (i.e. summed across all ranks, if distributed)
286 : * number of elements in a given block
287 : * @param[in] id subdomainID
288 : * return number of elements in block
289 : */
290 : unsigned int numElemsInSubdomain(const SubdomainID & id) const;
291 :
292 : /**
293 : * Whether the element is owned by this rank
294 : * @return whether element is owned by this rank
295 : */
296 : bool isLocalElem(const Elem * elem) const;
297 :
298 : /**
299 : * Get the global element ID from the local element ID
300 : * @param[in] id local element ID
301 : * @return global element ID
302 : */
303 15452357 : unsigned int globalElemID(const unsigned int & id) const { return _local_to_global_elem[id]; }
304 :
305 : /**
306 : * Set the cell temperature, and print helpful error message if a failure occurs; this sets
307 : * the temperature for the id + instance, which could be one of N cells filling the 'cell_info'
308 : * parent cell (which is what we actually use for error printing)
309 : * @param[in] id cell ID
310 : * @param[in] instance cell instance
311 : * @param[in] T temperature
312 : * @param[in] cell_info cell info for which we are setting interior temperature, for error printing
313 : */
314 : virtual void setCellTemperature(const int32_t & id, const int32_t & instance, const Real & T,
315 : const cellInfo & cell_info) const;
316 :
317 : /**
318 : * Set the cell density, and print helpful error message if a failure occurs
319 : * @param[in] density density
320 : * @param[in] cell_info cell info for which we are setting the density
321 : */
322 : virtual void setCellDensity(const Real & density, const cellInfo & cell_info) const;
323 :
324 : /**
325 : * Get a descriptive, formatted, string describing a cell
326 : * @param[in] cell_info cell index, instance pair
327 : * @param[in] brief whether to print a shorter string
328 : * @return descriptive string describing cell
329 : */
330 : virtual std::string printCell(const cellInfo & cell_info, const bool brief = false) const;
331 :
332 : /**
333 : * Get the fill of an OpenMC cell
334 : * @param[in] cell_info cell ID, instance
335 : * @param[out] fill_type fill type of the cell, one of MATERIAL, UNIVERSE, or LATTICE
336 : * @return indices of what is filling the cell
337 : */
338 : virtual std::vector<int32_t> cellFill(const cellInfo & cell_info, int & fill_type) const;
339 :
340 : /**
341 : * Whether the cell has a material fill (if so, then get the material index). Void counts
342 : * as a material, with a material index of -1.
343 : * @param[in] cell_info cell ID, instance
344 : * @param[out] material_index material index in the cell
345 : * @return whether the cell is filled by a material
346 : */
347 : bool materialFill(const cellInfo & cell_info, int32_t & material_index) const;
348 :
349 : /**
350 : * Calculate the number of unique OpenMC cells (each with a unique ID & instance)
351 : * @return number of unique OpenMC Cells in entire model
352 : */
353 : long unsigned int numCells() const;
354 :
355 : /**
356 : * Return all IDs of all Cardinal-mapped Tallies
357 : * @return all Cardinal-mapped Tally IDs
358 : */
359 : virtual std::vector<int32_t> getMappedTallyIDs() const = 0;
360 :
361 : /**
362 : * Whether or not the problem is tallying kinetics parameters.
363 : * @return if the problem is tallying kinetics parameters
364 : */
365 52 : bool computeKineticsParams() const { return _calc_kinetics_params; }
366 :
367 : /**
368 : * Get the tally responsible for accumulating the scores for \Lambda_{eff} and \beta_{eff}.
369 : * @return the global IFP tally
370 : */
371 : const openmc::Tally & getKineticsParamTally();
372 :
373 : protected:
374 : /// Find all userobjects which are changing OpenMC data structures
375 : void getOpenMCUserObjects();
376 :
377 : /// Ensure that the IDs of OpenMC objects in UserObjects don't clash
378 : void checkOpenMCUserObjectIDs() const;
379 :
380 : /// Ensure that any tally editors don't apply to Cardinal-mapped tallies
381 : void checkTallyEditorIDs() const;
382 :
383 : /// Execute all filter editor userobjects
384 : void executeFilterEditors();
385 :
386 : /// Execute all tally editor userobjects
387 : void executeTallyEditors();
388 :
389 : // execute tallly and filte editors
390 : void executeEditors();
391 :
392 : /// Set the nuclide densities for any materials being modified via MOOSE
393 : void sendNuclideDensitiesToOpenMC();
394 :
395 : /// Set the tally nuclides for any tallies being modified via MOOSE
396 : void sendTallyNuclidesToOpenMC();
397 :
398 : /**
399 : * Set an auxiliary elemental variable to a specified value
400 : * @param[in] var_num variable number
401 : * @param[in] elem_ids element IDs to set
402 : * @param[in] value value to set
403 : */
404 : void fillElementalAuxVariable(const unsigned int & var_num,
405 : const std::vector<unsigned int> & elem_ids,
406 : const Real & value);
407 :
408 : /**
409 : * Get name of source bank file to write
410 : * @param[out] file name
411 : */
412 40 : std::string sourceBankFileName() const
413 : {
414 80 : return openmc::settings::path_output + "initial_source_" +
415 80 : std::to_string(_fixed_point_iteration) + ".h5";
416 : }
417 :
418 : /// Whether to print diagnostic information about model setup and the transfers
419 : const bool & _verbose;
420 :
421 : /// Power by which to normalize the OpenMC results, for k-eigenvalue mode
422 : const Real * _power;
423 :
424 : /// Source strength by which to normalize the OpenMC results, for fixed source mode
425 : const Real * _source_strength;
426 :
427 : /**
428 : * Whether to take the starting fission source from iteration \f$n\f$ as the
429 : * fission source converged from iteration \f$n-1\f$. Setting this to true will
430 : * in most cases lead to improved convergence of the initial source as iterations
431 : * progress because you don't "start from scratch" each iteration and do the same
432 : * identical (within a random number seed) converging of the fission source.
433 : */
434 : bool _reuse_source;
435 :
436 : /// Whether a mesh scaling was specified by the user
437 : const bool _specified_scaling;
438 :
439 : /**
440 : * Multiplicative factor to apply to the mesh in the [Mesh] block in order to convert
441 : * whatever units that mesh is in into OpenMC's length scale of centimeters. For instance,
442 : * it is commonplace to develop fuel performance models with a length scale of meters.
443 : * Rather than onerously convert all MOOSE inputs to which OpenMC will be coupled to units
444 : * of centimeters, this parameter allows us to scale the mesh we couple to OpenMC.
445 : * Note that this does not actually scale the mesh itself, but simply multiplies the mesh
446 : * coordinates by this parameter when identifying the mapping between elements and cells.
447 : *
448 : * By default, this parameter is set to 1.0, meaning that OpenMC is coupled to another
449 : * MOOSE application with an input written in terms of centimeters. Set this parameter to 100
450 : * if the coupled MOOSE application is in units of meters, for instance.
451 : *
452 : * To summarize by example - if the MOOSE application uses units of meters, with a mesh
453 : * named mesh.exo, then the OpenMC-wrapped input file should also use that mesh (with
454 : * units of meters) in its [Mesh] block (or perhaps a coarser version of that mesh if
455 : * the resolution of coupling does not need to match - the units just have to be the same).
456 : * Then, you should set 'scaling = 100.0' so that the mapping is performed correctly.
457 : */
458 : const Real & _scaling;
459 :
460 : /// Whether to skip writing statepoints from OpenMC
461 : const bool & _skip_statepoint;
462 :
463 : /**
464 : * Fixed point iteration index used in relaxation; because we sometimes run OpenMC
465 : * in a pseudo-transient coupling with NekRS, we simply increment this by 1 each
466 : * time we call openmc::run(). This uses a zero indexing, so after the first iteration,
467 : * we have finished iteration 0, and so on.
468 : */
469 : int _fixed_point_iteration;
470 :
471 : /// Total number of particles simulated
472 : unsigned int _total_n_particles;
473 :
474 : /// Total number of unique OpenMC cell IDs + instances combinations
475 : long unsigned int _n_openmc_cells;
476 :
477 : /**
478 : * Number of digits to use to display the cell ID for diagnostic messages; this is
479 : * estimated conservatively based on the total number of cells, even though there
480 : * may be distributed cells such that the maximum cell ID is far smaller than the
481 : * total number of cells.
482 : */
483 : int _n_cell_digits;
484 :
485 : /// OpenMC run mode
486 : openmc::RunMode _run_mode;
487 :
488 : /// Userobjects for changing OpenMC material compositions
489 : std::vector<OpenMCNuclideDensities *> _nuclide_densities_uos;
490 :
491 : /// Userobjects for creating/changing OpenMC filters
492 : std::vector<OpenMCDomainFilterEditor *> _filter_editor_uos;
493 :
494 : /// Userobjects for creating/changing OpenMC tallies
495 : std::vector<OpenMCTallyEditor *> _tally_editor_uos;
496 :
497 : /// Mapping from local element indices to global element indices for this rank
498 : std::vector<unsigned int> _local_to_global_elem;
499 :
500 : /// Whether or not the problem contains mesh adaptivity.
501 : const bool _has_adaptivity;
502 :
503 : /**
504 : * A flag which is set to true if OpenMC should rerun after an adaptivity cycle.
505 : * This is a performance optimization to avoid rerunning OpenMC if adaptivity
506 : * didn't refine/coarsen the mesh on the last cycle.
507 : */
508 : bool _run_on_adaptivity_cycle;
509 :
510 : /// A flag to indicate if OpenMC is calculating kinetics parameters.
511 : const bool _calc_kinetics_params;
512 :
513 : /// Whether to reset the seed each time a new OpenMC calculation runs
514 : const bool & _reset_seed;
515 :
516 : /// The initial OpenMC seed
517 : const int64_t _initial_seed;
518 :
519 : /// The index of a global tally to accumulate the scores required for kinetics parameters.
520 : int _ifp_tally_index = -1;
521 :
522 : /// The global tally used to accumulate the scores required for kinetics parameters.
523 : openmc::Tally * _ifp_tally = nullptr;
524 :
525 : /// Conversion unit to transfer between kg/m3 and g/cm3
526 : static constexpr Real _density_conversion_factor{0.001};
527 :
528 : /// ID used by OpenMC to indicate that a material fill is VOID
529 : static constexpr int MATERIAL_VOID{-1};
530 : };
|
