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