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