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 : #include "OpenMCProblemBase.h"
22 : #include "SymmetryPointGenerator.h"
23 :
24 : /// Tally/filter includes.
25 : #include "TallyBase.h"
26 : #include "FilterBase.h"
27 :
28 : #ifdef ENABLE_DAGMC
29 : #include "MoabSkinner.h"
30 : #include "DagMC.hpp"
31 : #endif
32 :
33 : /// Forward declarations to avoid cyclic dependencies.
34 : class OpenMCVolumeCalculation;
35 :
36 : /**
37 : * Mapping of OpenMC to a collection of MOOSE elements, with temperature and/or
38 : * density feedback. The mappind is established authomatically by looping over
39 : * all the MOOSE elements and finding the OpenMC cell at each element's centroid.
40 : *
41 : * All feedback into OpenMC is performed via element averages. The
42 : * 'temperature_blocks' parameter is used to indicate which MOOSE blocks should
43 : * provide temperature feedback, while the 'density_blocks' parameter is used to
44 : * indicate which MOOSE blocks should provide density feedback. Tallies are
45 : * automatically added to OpenMC using either cell or mesh tallies.
46 : *
47 : * Each OpenMC cell shall not have ambiguous data transfers. That is, a cell
48 : * should map to a set of elements that are ALL/NONE providing temperature
49 : * feedback, ALL/NONE providing density feedback, ALL/NONE providing cell
50 : * tallies, and ALL/NONE being uncoupled altogether.
51 : *
52 : * TODO: If this is too restrictive in the future, we could implement some type
53 : * of weighted averaging process. Also, if a cell maps to a phase and an
54 : * unmapped region, perhaps we want to allow that.
55 : *
56 : * Other considerations you should be aware of:
57 : * - The density being transferred into OpenMC from MOOSE is in units of kg/m3; this
58 : * is the unit employed by the MOOSE fluid properties module.
59 : * - The temperature being transferred into OpenMC from MOOSE is in units of K; this
60 : * is the unit employed by the MOOSE fluid and solid properties modules.
61 : * - If your geometry is highly fine-scale (such as TRISOs), you might be able to get a
62 : * huge speedup in your runtime if you (i) build your OpenMC model by repeating the same
63 : * TRISO universe in each of your repeatable-units (e.g. pebbles, compacts, plates)
64 : * AND (ii) leverage the 'identical_cell_fills' option.
65 : */
66 : class OpenMCCellAverageProblem : public OpenMCProblemBase
67 : {
68 : public:
69 : OpenMCCellAverageProblem(const InputParameters & params);
70 : static InputParameters validParams();
71 :
72 : virtual void initialSetup() override;
73 : virtual void externalSolve() override;
74 : virtual void syncSolutions(ExternalProblem::Direction direction) override;
75 2679 : virtual bool converged(unsigned int) override { return true; }
76 :
77 : /**
78 : * Read a 2d vector of subdomain names, and check that there are no duplications
79 : * and that all provided values exist on the mesh.
80 : * @param[in] name string name for the 2d vector parameter
81 : * @param[out] names subdomain names
82 : * @param[out] flattened_ids flattened 1d vector of subdomain IDs
83 : */
84 : void read2DBlockParameters(const std::string name,
85 : std::vector<std::vector<SubdomainName>> & names,
86 : std::vector<SubdomainID> & flattened_ids);
87 :
88 : /**
89 : * Check that the specified blocks are in the mesh
90 : * @param[in] name name for throwing an error
91 : * @param[in] ids block IDs to check
92 : * @param[in] names block subdomain names for throwing an error
93 : */
94 : void checkBlocksInMesh(const std::string name,
95 : const std::vector<SubdomainID> & ids,
96 : const std::vector<SubdomainName> & names) const;
97 :
98 : /// Initialize the mapping of OpenMC to the MooseMesh and perform additional setup actions
99 : void setupProblem();
100 :
101 : /**
102 : * Add the tally variable(s) (to receive OpenMC tally values), temperature variable(s)
103 : * (to write into OpenMC cells), and density variable(s) (to write into OpenMC materials)
104 : */
105 : virtual void addExternalVariables() override;
106 :
107 : /**
108 : * Get the cell volume from a stochastic calculation
109 : * @param[in] cell_info cell index, instance pair
110 : * @return stochastically-computed OpenMC cell volume
111 : */
112 : virtual Real cellVolume(const cellInfo & cell_info) const;
113 :
114 : /**
115 : * Reference to stochastic volume calculation
116 : * @return reference to stochastic volume calculation
117 : */
118 609314 : virtual const OpenMCVolumeCalculation * volumeCalculation() const { return _volume_calc; }
119 :
120 : /**
121 : * Get the mapping of cells to MOOSE elements
122 : * @return mapping of cells to MOOSE elements
123 : */
124 12087 : virtual const std::map<cellInfo, std::vector<unsigned int>> & cellToElem() const
125 : {
126 12087 : return _cell_to_elem;
127 : }
128 :
129 : /**
130 : * Get the MOOSE subdomains associated with an OpenMC cell
131 : * @param info the cell info
132 : * @return MOOSE subdomains associated with an OpenMC cell
133 : */
134 14476 : virtual std::unordered_set<SubdomainID> getCellToElementSub(const cellInfo & info)
135 : {
136 14476 : return _cell_to_elem_subdomain.at(info);
137 : }
138 :
139 : /**
140 : * Whether transformations are applied to the [Mesh] points when mapping to OpenMC
141 : * @return whether transformations are applied
142 : */
143 2310221 : virtual bool hasPointTransformations() const { return _symmetry != nullptr; }
144 :
145 : /**
146 : * Get all the scores added to the tally
147 : * @return scores
148 : */
149 108 : virtual const std::vector<std::string> & getTallyScores() const { return _all_tally_scores; }
150 :
151 : /**
152 : * Get the number of tallies scoring a particular score.
153 : * @param[in] score the score to check
154 : * @return the number of tallies accumulating 'score'
155 : */
156 528 : unsigned int getNumScoringTallies(const std::string & score) const
157 : {
158 528 : return _score_count.count(score) > 0 ? _score_count.at(score) : 0;
159 : }
160 :
161 : /**
162 : * Check to see if this problem contains a specific tally score.
163 : * @param[in] score the tally score
164 : * @return whether this problem contains the tally score in a tally object
165 : */
166 356 : bool hasScore(const std::string & score)
167 : {
168 356 : return std::find(_all_tally_scores.begin(), _all_tally_scores.end(), score) !=
169 356 : _all_tally_scores.end();
170 : }
171 :
172 : /**
173 : * Get a tally object by its name.
174 : * @param[in] name the name of the TallyBase
175 : * @return a Cardinal wrapped tally
176 : */
177 : const TallyBase * getTally(const std::string & name);
178 :
179 : /**
180 : * Get the variable(s) associated with an OpenMC tally score.
181 : * @param[in] score the OpenMC score
182 : * @param[in] tally_name the name of the tally object to fetch score variables from
183 : * @param[in] tid the thread ID associated with the current MOOSE object
184 : * @param[in] output the output variable (relative error, standard deviation, etc.) to fetch
185 : * @param[in] skip_func_exp whether functional expansion filter bins should be skipped or not when
186 : * fetching variable values
187 : * @return a vector of variable values associated with score
188 : */
189 : std::vector<const MooseVariableFE<Real> *> getTallyScoreVariables(const std::string & score,
190 : const std::string & tally_name,
191 : THREAD_ID tid,
192 : const std::string & output = "",
193 : bool skip_func_exp = false);
194 :
195 : /**
196 : * Get the variable value(s) associated with an OpenMC tally score.
197 : * @param[in] score the OpenMC score
198 : * @param[in] tid the thread ID associated with the current MOOSE object
199 : * @param[in] tally_name the name of the tally object to fetch score variable values from
200 : * @param[in] output the output variable (relative error, standard deviation, etc.) to fetch
201 : * @param[in] skip_func_exp whether functional expansion filter bins should be skipped or not when
202 : * fetching variable values
203 : * @return a vector of variable values associated with score
204 : */
205 : std::vector<const VariableValue *> getTallyScoreVariableValues(const std::string & score,
206 : const std::string & tally_name,
207 : THREAD_ID tid,
208 : const std::string & output = "",
209 : bool skip_func_exp = false);
210 :
211 : /**
212 : * Get the variable value(s) associated with an OpenMC tally score.
213 : * @param[in] score the OpenMC score
214 : * @param[in] tally_name the name of the tally object to fetch score variable neighbor values from
215 : * @param[in] tid the thread ID associated with the current MOOSE object
216 : * @param[in] output the output variable (relative error, standard deviation, etc.) to fetch
217 : * @param[in] skip_func_exp whether functional expansion filter bins should be skipped or not when
218 : * fetching variable values
219 : * @return a vector of variable values associated with score
220 : */
221 : std::vector<const VariableValue *>
222 : getTallyScoreNeighborVariableValues(const std::string & score,
223 : const std::string & tally_name,
224 : THREAD_ID tid,
225 : const std::string & output = "",
226 : bool skip_func_exp = false);
227 :
228 : /**
229 : * Whether a tally contains a specified output or not.
230 : * @param[in] score the tally score to check
231 : * @param[in] output the additional output (unrelaxed standard deviation, relative error, or
232 : * tally)
233 : * @return whether an added tally has the output or not
234 : */
235 : bool hasOutput(const std::string & score, const std::string & output) const;
236 :
237 : /**
238 : * Apply transformations to point
239 : * @param[in] pt point
240 : * @return transformed point
241 : */
242 2261067 : virtual Point transformPoint(const Point & pt) const
243 : {
244 2261067 : return this->hasPointTransformations() ? _symmetry->transformPoint(pt) : pt;
245 : }
246 :
247 : /**
248 : * This class uses elem->volume() in order to normalize the tally values. However,
249 : * elem->volume() is expensive, so whenever MOOSE does integration, they set
250 : * _current_elem_volume to the volume as set by the sum of the quadrature weights.
251 : * The quadrature rule that MOOSE provides when you only have CONSTANT MONOMIALS is
252 : * insufficient for exactly integrating the element Jacobian mapping type (which
253 : * is FIRST LAGRANGE for a first order element), so you get an error relative to
254 : * the libmesh volume computation.
255 : *
256 : * So, we need to make sure that a minimum order quadrature rule is used
257 : * so that the total tally as computed by an
258 : * ElementIntegralVariablePostprocessor actually matches the specified total
259 : * (for low quadrature orders, there can be an error up to about 5% or so in total
260 : * power). This override simply forces the volume quadrature order to be 2 or higher
261 : * when using Gauss (default), monomial, or Gauss-Lobatto quadrature.
262 : *
263 : * For other quadrature rules, the approximations made in elem->volume() are never
264 : * going to match the volume integrations in MOOSE (unless the quadrature order is
265 : * very very high). For these orders, we print an error message informing the user
266 : * that they should switch to a different order.
267 : */
268 : virtual void createQRules(QuadratureType type,
269 : Order order,
270 : Order volume_order,
271 : Order face_order,
272 : SubdomainID block,
273 : bool allow_negative_weights = true) override;
274 :
275 : /**
276 : * Type definition for cells contained within a parent cell; the first value
277 : * is the cell index, while the second is the set of cell instances
278 : */
279 : typedef std::unordered_map<int32_t, std::vector<int32_t>> containedCells;
280 :
281 : /**
282 : * Get the cell index from the element ID; will return UNMAPPED for unmapped elements
283 : * @param[in] elem_id element ID
284 : * @return cell index
285 : */
286 : int32_t elemToCellIndex(const int & elem_id) const { return elemToCellInfo(elem_id).first; }
287 :
288 : /**
289 : * Get the cell ID from the element ID. Note that this function requires that the elem_id
290 : * maps to an OpenMC cell, or else an error will be raised from OpenMC in cellID.
291 : * @param[in] elem_id element ID
292 : * @return cell ID
293 : */
294 25794938 : int32_t elemToCellID(const int & elem_id) const { return cellID(elemToCellIndex(elem_id)); }
295 :
296 : /**
297 : * Get the cell instance from the element ID; will return UNMAPPED for unmapped elements
298 : * @param[in] elem_id element ID
299 : * @return cell instance
300 : */
301 : int32_t elemToCellInstance(const int & elem_id) const { return elemToCellInfo(elem_id).second; }
302 :
303 : /**
304 : * Get the cell index, instance pair from element ID; if the element doesn't map to an OpenMC
305 : * cell, the index and instance are both set to UNMAPPED
306 : * @param[in] elem_id element ID
307 : * @return cell index, instance pair
308 : */
309 37399002 : cellInfo elemToCellInfo(const int & elem_id) const { return _elem_to_cell[elem_id]; }
310 :
311 : /**
312 : * Get the cell material index based on index, instance pair. Note that this function requires
313 : * a valid instance, index pair for cellInfo - you cannot pass in an unmapped cell, i.e.
314 : * (UNMAPPED, UNMAPPED)
315 : * @param[in] cell_info cell index, instance pair
316 : * @return material index
317 : */
318 : int32_t cellToMaterialIndex(const cellInfo & cell_info) const
319 : {
320 2578688 : return _cell_to_material.at(cell_info);
321 : }
322 :
323 : /**
324 : * Get the fields coupled for each cell; because we require that each cell maps to a consistent
325 : * set, we simply look up the coupled fields of the first element that this cell maps to. Note
326 : * that this function requires a valid instance, index pair for cellInfo - you cannot pass in an
327 : * unmapped cell, i.e. (UNMAPPED, UNMAPPED)
328 : * @param[in] cell_info cell index, instance pair
329 : * @return coupling fields
330 : */
331 : coupling::CouplingFields cellFeedback(const cellInfo & cell_info) const;
332 :
333 : /**
334 : * Whether a cell has density feedback
335 : * @param[in] cell_info cell index,instance pair
336 : * @return if cell has density feedback
337 : */
338 5260047 : bool hasDensityFeedback(const cellInfo & cell_info) const
339 : {
340 : std::vector<coupling::CouplingFields> phase = {coupling::density,
341 5260047 : coupling::density_and_temperature};
342 5260047 : return std::find(phase.begin(), phase.end(), cellFeedback(cell_info)) != phase.end();
343 5260047 : }
344 :
345 : /**
346 : * Whether a cell has temperature feedback
347 : * @param[in] cell_info cell index,instance pair
348 : * @return if cell has temperature feedback
349 : */
350 4093276 : bool hasTemperatureFeedback(const cellInfo & cell_info) const
351 : {
352 : std::vector<coupling::CouplingFields> phase = {coupling::temperature,
353 4093276 : coupling::density_and_temperature};
354 4093276 : return std::find(phase.begin(), phase.end(), cellFeedback(cell_info)) != phase.end();
355 4093276 : }
356 :
357 : /**
358 : * Checks if the [Problem/Filters] block contains a specific filter.
359 : * @param[in] filter_name the MOOSE object name of the filter
360 : * @return whether the problem contains the specified filter
361 : */
362 : bool hasFilter(const std::string & filter_name) const { return _filters.count(filter_name) > 0; }
363 :
364 : /**
365 : * Get a filter added by the [Problem/Filters] block by it's MOOSE object name.
366 : * @param[in] filter_name the MOOSE object name of the filter
367 : * @return the filter object
368 : */
369 : std::shared_ptr<FilterBase> & getFilter(const std::string & filter_name)
370 : {
371 1054 : return _filters.at(filter_name);
372 : }
373 :
374 : /**
375 : * Get the local tally
376 : * @return local tally
377 : */
378 : const std::vector<std::shared_ptr<TallyBase>> & getLocalTallies() const { return _local_tallies; }
379 :
380 : /**
381 : * Get the temperature of a cell; for cells not filled with materials, this will return
382 : * the temperature of the first material-type cell
383 : * @param[in] cell_info cell index, instance pair
384 : */
385 : double cellTemperature(const cellInfo & cell_info) const;
386 :
387 : /**
388 : * Get the volume that each OpenMC cell mapped to
389 : * @param[in] cell_info cell index, instance pair
390 : */
391 : double cellMappedVolume(const cellInfo & cell_info) const;
392 :
393 : /// Reconstruct the DAGMC geometry after skinning
394 : void reloadDAGMC();
395 :
396 : /**
397 : * Add a Filter object using the filter system.
398 : * @param[in] type the new tally type
399 : * @param[in] name the name of the new tally
400 : * @param[in] moose_object_pars the input parameters of the new tally
401 : */
402 : void addFilter(const std::string & type,
403 : const std::string & name,
404 : InputParameters & moose_object_pars);
405 :
406 : /**
407 : * Add a Tally object using the tally system.
408 : * @param[in] type the new tally type
409 : * @param[in] name the name of the new tally
410 : * @param[in] moose_object_pars the input parameters of the new tally
411 : * @return a shared pointer to the TallyBase
412 : */
413 : std::shared_ptr<TallyBase>
414 : addTally(const std::string & type, const std::string & name, InputParameters & moose_object_pars);
415 :
416 : /**
417 : * Multiplier on the normalized tally results; for fixed source runs,
418 : * we multiply the tally (which has units of eV/source)
419 : * by the source strength and the eV to joule conversion, while for k-eigenvalue runs, we
420 : * multiply the normalized tally (which is unitless and has an integral
421 : * value of 1.0) by the power.
422 : * @param[in] score_name name of the score
423 : * @param[in] local_mean_tally the mean tally associated with score_name
424 : */
425 : Real tallyMultiplier(const std::string & score_name, const Real & local_mean_tally) const;
426 :
427 : /**
428 : * Get the reference density of an element when running in multi-group mode.
429 : * @param[in] elem the element
430 : * @return the reference density (kg/m3) or unity (if not running in multi-group mode)
431 : */
432 : const Real getReferenceDensity(const Elem * elem) const;
433 :
434 : /**
435 : * Check whether a vector extracted with getParam is empty
436 : * @param[in] vector vector
437 : * @param[in] name name to use for printing error if empty
438 : */
439 : template <typename T>
440 4672 : void checkEmptyVector(const std::vector<T> & vector, const std::string & name) const
441 : {
442 4566 : if (vector.empty())
443 8 : mooseError(name + " cannot be empty!");
444 4664 : }
445 :
446 4054 : int fixedPointIteration() const { return _fixed_point_iteration; }
447 :
448 : /**
449 : * Checks if the problem uses adaptivity or not.
450 : * @return if the problem uses adaptivity.
451 : */
452 2839 : bool hasAdaptivity() const { return _has_adaptivity; }
453 :
454 : /**
455 : * Checks if the problem is using a MoabSkinner or not.
456 : * @return if the problem uses skinning
457 : */
458 157 : bool hasSkinner() const { return _using_skinner; }
459 :
460 : /// Constant flag to indicate that a cell/element was unmapped
461 : static constexpr int32_t UNMAPPED{-1};
462 :
463 : /// Spatial dimension of the Monte Carlo problem
464 : static constexpr int DIMENSION{3};
465 :
466 : /// Get a modifyable non-const reference to the Moose mesh
467 : virtual MooseMesh & getMooseMesh();
468 :
469 : /// Get a modifyable const reference to the Moose mesh
470 : virtual const MooseMesh & getMooseMesh() const;
471 :
472 : /**
473 : * Whether a moving mesh is used
474 : * @return whether the [Mesh] is moving
475 : */
476 : const bool & useDisplaced() const { return _use_displaced; }
477 :
478 : protected:
479 : /**
480 : * A function to re-initialize coupling and apply feedback to the OpenMC problem.
481 : * Applied before OpenMC is executed in either: i) normal Picard iterations, or
482 : * ii) criticality searches. This function performs these operations in the
483 : * following order:
484 : * 1. Updates the OpenMC geometry using the skinner;
485 : * 2. Resets tallies and the cell->element maps to take into account
486 : * mesh/geometry changes;
487 : * 3. Updates the nuclide composition of OpenMC materials;
488 : * 4. Sets cell temperatures and densities;
489 : * 5. Exports OpenMC properties;
490 : * 6. Reinitializes multi-group cross sections to take into account changing
491 : * temperatures and densities.
492 : */
493 : void reinitCouplingAndApplyFeedback();
494 :
495 : /**
496 : * Implement critSearchStep() to re-generate the cell-to-element (and dual)
497 : * mapping. This sends new temperatures and densities to OpenMC from the
498 : * re-mapped elements to ensure the state remains critical under changes to the
499 : * model with feedbacks.
500 : */
501 : virtual void critSearchStep() override;
502 :
503 : /**
504 : * Get the cell level in OpenMC to use for coupling
505 : * @param[in] c point
506 : * @return cell level
507 : */
508 : unsigned int getCellLevel(const Point & c) const;
509 :
510 : /**
511 : * Read the names of the MOOSE variables used for sending feedback into OpenMC
512 : * @param[in] param feedback term to read
513 : * @param[in] default_name default name to use for MOOSE variables holding this field
514 : * @param[out] vars_to_specified_blocks map from MOOSE variable names to the blocks on which they
515 : * are defined
516 : * @param[out] specified_blocks user-specified blocks for feedback
517 : */
518 : void
519 : readBlockVariables(const std::string & param,
520 : const std::string & default_name,
521 : std::map<std::string, std::vector<SubdomainName>> & vars_to_specified_blocks,
522 : std::vector<SubdomainID> & specified_blocks);
523 :
524 : /**
525 : * Whether this cell has an identical fill
526 : * @param[in] cell_info cell index, instance pair
527 : * @return whether this cell has an identical fill
528 : */
529 : bool cellHasIdenticalFill(const cellInfo & cell_info) const;
530 :
531 : /**
532 : * When using the 'identical_cell_fills' feature, this is used to determine the
533 : * contained material cells in each parent cell by applying a uniform shift
534 : * @param[in] cell_info cell index, instance pair
535 : * @return material cells contained within the given cell
536 : */
537 : containedCells shiftCellInstances(const cellInfo & cell_info) const;
538 :
539 : /**
540 : * Whether this cell overlaps with ANY value in the given subdomain set
541 : * @param[in] cell_info cell index, instance pair
542 : * @param[in] id subdomain IDs
543 : * @return whether the cell overlaps with the subdomain
544 : */
545 : bool cellMapsToSubdomain(const cellInfo & cell_info,
546 : const std::unordered_set<SubdomainID> & id) const;
547 :
548 : /**
549 : * Get the first material cell contained in the given cell
550 : * @param[in] cell_info cell index, instance pair
551 : * @return material cell index, instance pair
552 : */
553 : cellInfo firstContainedMaterialCell(const cellInfo & cell_info) const;
554 :
555 : /**
556 : * Get all of the material cells contained within this cell
557 : * @param[in] cell_info cell index, instance pair
558 : * @return all material cells contained in the given cell
559 : */
560 : containedCells containedMaterialCells(const cellInfo & cell_info) const;
561 :
562 : /**
563 : * Delete the OpenMC DAGMC geometry and re-generate the CSG geometry data structures in-place.
564 : */
565 : void updateOpenMCGeometry();
566 :
567 : /**
568 : * Get a list of each material in the problem, sorted by subdomain. This function also checks
569 : * that there is just one OpenMC material in each subdomain, necessary for the DAGMC skinning.
570 : * @return material in each subdomain
571 : */
572 : std::vector<std::string> getMaterialInEachSubdomain() const;
573 :
574 : /**
575 : * Apply transformations and scale point from MOOSE into the OpenMC domain
576 : * @param[in] pt point
577 : * @return transformed point
578 : */
579 : Point transformPointToOpenMC(const Point & pt) const;
580 :
581 : /**
582 : * For geometries with fine-scale details (e.g. TRISO), Cardinal's default settings can
583 : * take a very long time to initialize the problem (but we can't change those defaults
584 : * because they are not 100% applicable all the time). So, we print out a message to
585 : * the user to point them in the right direction if their initialization is taking a
586 : * long time.
587 : * @param[in] start time to use for evaluating whether we've exceeded our limit for printing the
588 : * message
589 : */
590 : void
591 : printTrisoHelp(const std::chrono::time_point<std::chrono::high_resolution_clock> & start) const;
592 :
593 : /**
594 : * Print to the console the names of the auxvariables used for I/O with OpenMC.
595 : * We only print these tables once, upon initialization, because this data does
596 : * not change if re-initializing the spatial mapping for moving-mesh problems,
597 : * adaptive refinement, skinning, etc.
598 : */
599 : void printAuxVariableIO();
600 :
601 : /**
602 : * Get all the material indices within the set of cells
603 : * @param[in] contained_cells set of cells
604 : * @return contained materials
605 : */
606 : std::vector<int32_t> materialsInCells(const containedCells & contained_cells) const;
607 :
608 : /// Loop over the mapped cells, and build a map between subdomains to OpenMC materials
609 : void subdomainsToMaterials();
610 :
611 : /**
612 : * Get a set of all subdomains that have at least 1 element coupled to an OpenMC cell
613 : * @return subdomains with at least 1 element coupled to OpenMC
614 : */
615 : std::set<SubdomainID> coupledSubdomains() const;
616 :
617 : /**
618 : * Gather a vector of values to be summed for each cell
619 : * @param[in] local local values to be summed for the cells
620 : * @param[out] global global mapping of the summed values to the cells
621 : */
622 : template <typename T>
623 : void gatherCellSum(std::vector<T> & local, std::map<cellInfo, T> & global) const;
624 :
625 : /**
626 : * Gather a vector of values to be pushed back to for each cell
627 : * @param[in] local local values to be pushed back for the cells
628 : * @param[in] n_local number of local values contributed to each cell
629 : * @param[out] global global mapping of the pushed back values to the cells
630 : */
631 : template <typename T>
632 : void gatherCellVector(std::vector<T> & local, std::vector<unsigned int> & n_local, std::map<cellInfo, std::vector<T>> & global);
633 :
634 : /**
635 : * Get the feedback which this element provides to OpenMC
636 : * @param[in] elem
637 : * @return coupling phase
638 : */
639 : coupling::CouplingFields elemFeedback(const Elem * elem) const;
640 :
641 : /**
642 : * Read the parameters needed for triggers
643 : * @param[in] params input parameters
644 : */
645 : void getTallyTriggerParameters(const InputParameters & params);
646 :
647 : /**
648 : * Read the block parameters based on user settings
649 : * @param[in] name name of input parameter representing a vector of subdomain names
650 : * @param[in] blocks list of block ids to write
651 : */
652 : void readBlockParameters(const std::string name, std::unordered_set<SubdomainID> & blocks);
653 :
654 : /**
655 : * Cache the material cells contained within each coupling cell;
656 : * depending on user settings, this may attempt to take shortcuts
657 : * by assuming each cell has the same fills
658 : */
659 : void cacheContainedCells();
660 :
661 : /**
662 : * Fill the cached contained cells data structure for a given cell
663 : * @param[in] cell_info cell index, instance pair
664 : * @param[in] hint location hint used to accelerate the search
665 : * @param[out] map contained cell map
666 : */
667 : void setContainedCells(const cellInfo & cell_info,
668 : const Point & hint,
669 : std::map<cellInfo, containedCells> & map);
670 :
671 : /**
672 : * Check that the structure of the contained material cells for two cell matches;
673 : * i.e. this checks that the keys are the same and that the *number* of instances
674 : * of each filling material cell match.
675 : * @param[in] cell_info cell index, instance pair
676 : * @param[in] reference map we want to check against
677 : * @param[in] compare map we want to check
678 : */
679 : void checkContainedCellsStructure(const cellInfo & cell_info,
680 : containedCells & reference,
681 : containedCells & compare) const;
682 :
683 : /**
684 : * Set a minimum order for a volume quadrature rule
685 : * @param[in] volume_order order of the volume quadrature rule
686 : * @param[in] type string type of quadrature rule for printing a console message
687 : */
688 : void setMinimumVolumeQRules(Order & volume_order, const std::string & type);
689 :
690 : /// For keeping the output neat when using verbose
691 : std::string printNewline() const
692 : {
693 : if (_verbose)
694 : return "\n";
695 : else
696 : return "";
697 : }
698 :
699 : /// Loop over the elements in the MOOSE mesh and store the type of feedback applied by each.
700 : void storeElementPhase();
701 :
702 : /**
703 : * Loop over all the OpenMC cells and count the number of MOOSE elements to which the cell
704 : * is mapped based on phase.
705 : */
706 : void getCellMappedPhase();
707 :
708 : /// This function is used to ensure that each OpenMC cell only maps to a single phase
709 : void checkCellMappedPhase();
710 :
711 : /// Loop over all the OpenMC cells and get the element subdomain IDs that map to each cell
712 : void getCellMappedSubdomains();
713 :
714 : /**
715 : * Loop over all the OpenMC cells and compute the volume of the MOOSE elements that each
716 : * cell maps to
717 : */
718 : void computeCellMappedVolumes();
719 :
720 : /// Set up the mapping from MOOSE elements to OpenMC cells
721 : void initializeElementToCellMapping();
722 :
723 : /// Populate maps of MOOSE elements to OpenMC cells
724 : void mapElemsToCells();
725 :
726 : /**
727 : * A function which validates local tallies. This is done to ensure that at least one of the
728 : * tallies contains a heating score when running in eigenvalue mode. This must be done outside
729 : * of the constructor as tallies are added from an external system.
730 : */
731 : void validateLocalTallies();
732 :
733 : /// Add OpenMC tallies to facilitate the coupling
734 : void initializeTallies();
735 :
736 : /**
737 : * Reset any tallies previously added by Cardinal, by deleting them from OpenMC.
738 : * Also delete any mesh filters and meshes added to OpenMC for mesh filters.
739 : */
740 : void resetTallies();
741 :
742 : /// Find the material filling each cell which receives density feedback
743 : void getMaterialFills();
744 :
745 : /**
746 : * Get one point inside each cell, for accelerating the particle search routine.
747 : * This function will get the centroid of the first global element in the lowest
748 : * rank in the cell.
749 : */
750 : void getPointInCell();
751 :
752 : /**
753 : * Compute the product of volume with a field across ranks and sum into a global map
754 : * @param[in] var_num variable to weight with volume, mapped by subdomain ID
755 : * @param[in] phase phases to compute the operation for
756 : * @param[in] scaling a scaling factor to apply, mapped by subdomain ID
757 : * @return volume-weighted field for each cell, in a global sense
758 : */
759 : std::map<cellInfo, Real> computeVolumeWeightedCellInput(
760 : const std::map<SubdomainID, std::pair<unsigned int, std::string>> & var_num,
761 : const std::vector<coupling::CouplingFields> * phase = nullptr,
762 : const std::map<SubdomainID, Real> * scaling = nullptr) const;
763 :
764 : /**
765 : * Send temperature from MOOSE to OpenMC by computing a volume average
766 : * and applying a single temperature per OpenMC cell
767 : */
768 : void sendTemperatureToOpenMC() const;
769 :
770 : /**
771 : * Send density from MOOSE to OpenMC by computing a volume average
772 : * and applying a single density per OpenMC cell.
773 : */
774 : void sendDensityToOpenMC() const;
775 :
776 : /**
777 : * Check if a mapped location is in the outer universe of a lattice
778 : * @param[in] level lattice level
779 : * @return whether the location is in the outer universe
780 : */
781 : void latticeOuterCheck(const Point & c, int level) const;
782 :
783 : /**
784 : * Report an error for a mapped location in an outer universe of a lattice
785 : * @param[in] c Mapped location
786 : * @param[in] level level of the mapped cell
787 : */
788 : void latticeOuterError(const Point & c, int level) const;
789 :
790 : /**
791 : * Find the OpenMC cell at a given point in space
792 : * @param[in] point point
793 : * @return whether OpenMC reported an error
794 : */
795 : bool findCell(const Point & point);
796 :
797 : /**
798 : * Checks that the contained material cells exactly match between a reference obtained
799 : * by calling openmc::Cell::get_contained_cells for each cell and a shortcut
800 : * approach that assumes all identical cells (which aren't simply just material fills)
801 : * has exactly the same contained material cells.
802 : * @param[in] reference reference map to compare against
803 : * @param[in] compare shortcut map to compare
804 : */
805 : void compareContainedCells(std::map<cellInfo, containedCells> & reference,
806 : std::map<cellInfo, containedCells> & compare) const;
807 :
808 : /**
809 : * Return all IDs of all Cardinal-mapped Tallies
810 : * @return all Cardinal-mapped Tally IDs
811 : */
812 : virtual std::vector<int32_t> getMappedTallyIDs() const override;
813 :
814 : /// A reference to the serialized auxvariable solution.
815 : NumericVector<Number> & _serialized_solution;
816 :
817 : /**
818 : * Whether to automatically compute the mapping of OpenMC cell IDs and
819 : * instances to the [Mesh].
820 : */
821 : const bool & _output_cell_mapping;
822 :
823 : /**
824 : * Where to get the initial OpenMC temperatures and densities from;
825 : * can be either hdf5 (from a properties.h5 file), xml (whatever is already
826 : * set in the XML files), or moose (meaning whatever ICs are set on the 'temperature_variables'
827 : * and 'density_variables'
828 : */
829 : const coupling::OpenMCInitialCondition _initial_condition;
830 :
831 : /// Type of relaxation to apply to the OpenMC tallies
832 : const relaxation::RelaxationEnum _relaxation;
833 :
834 : /**
835 : * Type of trigger to apply to k eigenvalue to indicate when
836 : * the simulation is complete. These can be used to on-the-fly adjust the number
837 : * of active batches in order to reach some desired criteria (which is specified
838 : * by this parameter).
839 : */
840 : const trigger::TallyTriggerTypeEnum _k_trigger;
841 :
842 : /**
843 : * Coordinate level in the OpenMC domain to use for mapping cells to mesh.
844 : * When using 'lowest_cell_level', this parameter indicates that the lowest
845 : * cell level is used, up until _cell_level.
846 : */
847 : unsigned int _cell_level;
848 :
849 : /**
850 : * Whether OpenMC properties (temperature and density) should be exported
851 : * after being updated in syncSolutions.
852 : */
853 : const bool & _export_properties;
854 :
855 : /// Whether or not the problem uses a skinner to regenerate the OpenMC geometry.
856 : const bool _using_skinner;
857 :
858 : /**
859 : * When the mesh changes during the simulation (either from adaptive mesh refinement
860 : * or deformation), the mapping from OpenMC cells to the [Mesh] must be re-established
861 : * after each OpenMC run.
862 : */
863 : bool _need_to_reinit_coupling;
864 :
865 : /**
866 : * If known a priori by the user, whether the tally cells (which are not simply material
867 : * fills) have EXACTLY the same contained material cells. This is a big optimization for
868 : * TRISO problems in setting up homogenized temperature/density feedback to OpenMC.
869 : *
870 : * The concept can best be explained with a pebble bed reactor.
871 : * If every pebble is filled with an identical TRISO universe, then the material fills
872 : * in each pebble are identical to one another except for a constant offset. This idea
873 : * can be used to then skip all but the first two openmc::Cell::get_contained_cells
874 : * calls (which are required in order to figure out the pattern by which pebble N is
875 : * incremented relative to pebble 1).
876 : *
877 : * When using this parameter, we HIGHLY recommend setting 'check_identical_cell_fills =
878 : * true' the first time you run your model. This will figure out the material cell fills using a
879 : * method that calls openmc::Cell::get_contained_cells for every tally cell, i.e. without assuming
880 : * anything about repeated structure in your OpenMC model. Setting 'identical_cell_fills'
881 : * without also setting 'check_identical_cell_fills = true' may result in SILENT
882 : * errors!!! So it is essential to be sure you've removed any error sources before you turn the
883 : * error check off to actually leverage the speedup.
884 : *
885 : * Note: for any tally cells that are just filled with a material, we use the approach
886 : * where openmc::Cell::get_contained_cells is called in full.
887 : *
888 : * This optimization will not work (and 'check_identical_cell_fills = true' *will*
889 : * catch these) for:
890 : * - any situation where tallied, non-material-fill pebbles have different fills
891 : * (such as if you have different TRISO lattices in each pebble)
892 : * - any situation where there is a "gap" in the incrementing of the material fill
893 : * instances (such as if pebble 89 does not map to 'tally_blocks', then the instance
894 : * shift for pebble 90 relative to pebble 1 is 89, when it should have been 90).
895 : */
896 : const bool _has_identical_cell_fills;
897 :
898 : /**
899 : * Whether we should rigorously check that each tally cell has identical fills;
900 : * this is SLOW for large TRISO problems, but is essential to ensure the accuracy of
901 : * 'identical_cell_fills'. Please set 'check_identical_cell_fills' to 'true' at least
902 : * once before running production cases to be sure the optimization can be applied.
903 : */
904 : const bool & _check_identical_cell_fills;
905 :
906 : /**
907 : * Whether it can be assumed that all of the tallies (both those set by the user
908 : * in the XML file, as well as those created automatically by Cardinal) are
909 : * spatially separate. This means that once a particle scores to one tally bin, it wouldn't
910 : * score to ANY other tally bins. This can dramatically increase tracking rates
911 : * for problems with many tallies.
912 : */
913 : const bool & _assume_separate_tallies;
914 :
915 : /**
916 : * Whether to map density according to each individual OpenMC cell (in which case an
917 : * error is thrown if you don't have a unique material in each cell) or by material.
918 : */
919 : bool _map_density_by_cell;
920 :
921 : /**
922 : * Whether the problem has density feedback blocks specified; note that this is NOT necessarily
923 : * indicative that the mapping was successful in finding any cells corresponding to those blocks
924 : */
925 : const bool _specified_density_feedback;
926 :
927 : /**
928 : * Whether the problem has temperature feedback blocks specified; note that this is NOT
929 : * necessarily indicative that the mapping was successful in finding any cells corresponding to
930 : * those blocks
931 : */
932 : const bool _specified_temperature_feedback;
933 :
934 : /// Whether any cell tallies exist.
935 : bool _has_cell_tallies = false;
936 :
937 : /// Whether any spatial mapping from OpenMC's cells to the mesh is needed
938 : bool _needs_to_map_cells;
939 :
940 : /**
941 : * A map of the filter objects created by the [Problem/Filters] block. The key for each filter is
942 : * it's corresponding MOOSE name to allow tallies to look up filters.
943 : */
944 : std::map<std::string, std::shared_ptr<FilterBase>> _filters;
945 :
946 : /// A vector of the tally objects created by the [Problem/Tallies] block.
947 : std::vector<std::shared_ptr<TallyBase>> _local_tallies;
948 :
949 : /// A list of all of the scores contained by the local tallies added in the [Tallies] block.
950 : std::vector<std::string> _all_tally_scores;
951 :
952 : /// Number of tallies scoring a particular score.
953 : std::map<std::string, unsigned int> _score_count;
954 :
955 : /// A vector of auxvariable ids added by the [Tallies] block.
956 : std::vector<std::vector<unsigned int>> _tally_var_ids;
957 :
958 : /// A vector of external (output-based) auxvariable ids added by the [Tallies] block.
959 : std::vector<std::vector<std::vector<unsigned int>>> _tally_ext_var_ids;
960 :
961 : /// Blocks in MOOSE mesh that provide density feedback
962 : std::vector<SubdomainID> _density_blocks;
963 :
964 : /// Blocks in MOOSE mesh that provide temperature feedback
965 : std::vector<SubdomainID> _temp_blocks;
966 :
967 : /// Blocks for which the cell fills are identical
968 : std::unordered_set<SubdomainID> _identical_cell_fill_blocks;
969 :
970 : /// Mapping of MOOSE elements to the OpenMC cell they map to (if any)
971 : std::vector<cellInfo> _elem_to_cell{};
972 :
973 : /// Phase of each cell
974 : std::map<cellInfo, coupling::CouplingFields> _cell_phase;
975 :
976 : /// Number of elements in the MOOSE mesh that exclusively provide density feedback
977 : int _n_moose_density_elems;
978 :
979 : /// Number of elements in the MOOSE mesh that exclusively provide temperature feedback
980 : int _n_moose_temp_elems;
981 :
982 : /// Number of elements in the MOOSE mesh which provide temperature+density feedback
983 : int _n_moose_temp_density_elems;
984 :
985 : /// Number of no-coupling elements in the MOOSE mesh
986 : int _n_moose_none_elems;
987 :
988 : /**
989 : * Number of MOOSE elements that exclusively provide temperature feedback,
990 : * and which successfully mapped to OpenMC cells
991 : */
992 : int _n_mapped_temp_elems;
993 :
994 : /**
995 : * Number of MOOSE elements that exclusively provide density feedback,
996 : * and which successfully mapped to OpenMC cells
997 : */
998 : int _n_mapped_density_elems;
999 :
1000 : /**
1001 : * Number of MOOSE elements that provide temperature+density feedback,
1002 : * and which successfully mapped to OpenMC cells
1003 : */
1004 : int _n_mapped_temp_density_elems;
1005 :
1006 : /// Number of no-coupling elements mapped to OpenMC cells
1007 : int _n_mapped_none_elems;
1008 :
1009 : /// Total volume of uncoupled MOOSE mesh elements
1010 : Real _uncoupled_volume;
1011 :
1012 : /// Whether non-material cells are mapped
1013 : bool _material_cells_only{true};
1014 :
1015 : /// Mapping of OpenMC cell indices to a vector of MOOSE element IDs
1016 : std::map<cellInfo, std::vector<unsigned int>> _cell_to_elem;
1017 :
1018 : /// Mapping of OpenMC cell indices to a vector of MOOSE element IDs, on each local rank
1019 : std::map<cellInfo, std::vector<unsigned int>> _local_cell_to_elem;
1020 :
1021 : /// Mapping of OpenMC cell indices to the subdomain IDs each maps to
1022 : std::map<cellInfo, std::unordered_set<SubdomainID>> _cell_to_elem_subdomain;
1023 :
1024 : /// Mapping of elem subdomains to materials
1025 : std::map<SubdomainID, std::set<int32_t>> _subdomain_to_material;
1026 :
1027 : /**
1028 : * A point inside the cell, taken simply as the centroid of the first global
1029 : * element inside the cell. This is stored to accelerate the particle search.
1030 : */
1031 : std::map<cellInfo, Point> _cell_to_point;
1032 :
1033 : /**
1034 : * Volume associated with the mapped element space for each OpenMC cell; the unit
1035 : * for this volume is whatever is used in the [Mesh] block
1036 : */
1037 : std::map<cellInfo, Real> _cell_to_elem_volume;
1038 :
1039 : /**
1040 : * Volume associated with the actual OpenMC cell, computed by an optional
1041 : * OpenMCVolumeCalculation user object
1042 : */
1043 : std::map<cellInfo, Real> _cell_volume;
1044 :
1045 : /**
1046 : * Material filling each cell to receive density feedback. We enforce that these
1047 : * cells are filled with a material (cannot be filled with a lattice or universe).
1048 : */
1049 : std::map<cellInfo, int32_t> _cell_to_material;
1050 :
1051 : /**
1052 : * Material-type cells contained within a cell; this is only populated if a cell
1053 : * is NOT indicated as having an identical fill
1054 : */
1055 : std::map<cellInfo, containedCells> _cell_to_contained_material_cells;
1056 :
1057 : /// Number of material-type cells contained within a cell
1058 : std::map<cellInfo, int32_t> _cell_to_n_contained;
1059 :
1060 : /// Whether the present transfer is the first transfer
1061 : static bool _first_transfer;
1062 :
1063 : /// Whether the diagnostic tables on initialization have already been printed
1064 : static bool _printed_initial;
1065 :
1066 : /// Whether a warning has already been printed about very long setup times (for TRISOs)
1067 : static bool _printed_triso_warning;
1068 :
1069 : /// Dummy particle to reduce number of allocations of particles for cell lookup routines
1070 : openmc::Particle _particle;
1071 :
1072 : /// Number of particles simulated in the first iteration in Dufek-Gudowski relaxation
1073 : unsigned int _n_particles_1;
1074 :
1075 : /// Mapping from temperature variable name to the subdomains on which to read it from
1076 : std::map<std::string, std::vector<SubdomainName>> _temp_vars_to_blocks;
1077 :
1078 : /// Mapping from density variable name to the subdomains on which to read it from
1079 : std::map<std::string, std::vector<SubdomainName>> _density_vars_to_blocks;
1080 :
1081 : /// Optional volume calculation for cells which map to MOOSE
1082 : OpenMCVolumeCalculation * _volume_calc;
1083 :
1084 : /// Userobject that maps from a partial-symmetry OpenMC model to a whole-domain [Mesh]
1085 : const SymmetryPointGenerator * _symmetry;
1086 :
1087 : /// Number of temperature-only feedback elements in each mapped OpenMC cell (global)
1088 : std::map<cellInfo, int> _n_temp;
1089 :
1090 : /// Number of density-only feedback elements in each mapped OpenMC cell (global)
1091 : std::map<cellInfo, int> _n_rho;
1092 :
1093 : /// Number of temperature+density feedback elements in each mapped OpenMC cell (global)
1094 : std::map<cellInfo, int> _n_temp_rho;
1095 :
1096 : /// Number of none elements in each mapped OpenMC cell (global)
1097 : std::map<cellInfo, int> _n_none;
1098 :
1099 : /// The tally to be used for normalizing all other tallies when running an eigenvalue calculation.
1100 : std::shared_ptr<TallyBase> _source_rate_norm_tally;
1101 :
1102 : /// The score index into "_source_rate_norm_tally".
1103 : int _source_rate_score = -1;
1104 :
1105 : #ifdef ENABLE_DAGMC
1106 : /// Optional skinner to re-generate the OpenMC geometry on-the-fly for DAGMC models
1107 : MoabSkinner * _skinner = nullptr;
1108 :
1109 : /// Pointer to DAGMC
1110 : std::shared_ptr<moab::DagMC> _dagmc = nullptr;
1111 : #endif
1112 :
1113 : /// Total number of unique OpenMC cell IDs + instances combinations
1114 : long unsigned int _n_openmc_cells;
1115 :
1116 : /// ID of the OpenMC universe corresponding to the DAGMC universe
1117 : int32_t _dagmc_universe_id;
1118 :
1119 : /// Whether the DAGMC universe is the root universe or not.
1120 : bool _dagmc_root_universe = true;
1121 :
1122 : /// ID of the OpenMC cell corresponding to the cell which uses the DAGMC universe as a fill.
1123 : int32_t _cell_using_dagmc_universe_id;
1124 :
1125 : /// The number of OpenMC surfaces before skinning occurs. This is required to properly reinitialize
1126 : /// the CSG geometry contained in the OpenMC model.
1127 : const int32_t _initial_num_openmc_surfaces;
1128 :
1129 : /// Conversion rate from eV to Joule
1130 : static constexpr Real EV_TO_JOULE = 1.6022e-19;
1131 :
1132 : /// Tolerance for setting zero tally
1133 : static constexpr Real ZERO_TALLY_THRESHOLD = 1e-12;
1134 :
1135 : private:
1136 : /**
1137 : * Update the number of particles according to the Dufek-Gudowski relaxation scheme
1138 : */
1139 : void dufekGudowskiParticleUpdate();
1140 :
1141 : /// Flattened cell IDs collected after parallel communication
1142 : std::vector<int32_t> _flattened_ids;
1143 :
1144 : /// Flattened cell instancess collected after parallel communication
1145 : std::vector<int32_t> _flattened_instances;
1146 :
1147 : /// Offsets for each cell instance in an identically-repeated universe
1148 : containedCells _instance_offsets;
1149 :
1150 : /// Offset for each cell relative to the first identical-fill cell
1151 : std::map<cellInfo, int32_t> _n_offset;
1152 :
1153 : /// First identical-fill cell
1154 : cellInfo _first_identical_cell;
1155 :
1156 : /// Materials in the first identical-fill cell
1157 : std::vector<int32_t> _first_identical_cell_materials;
1158 :
1159 : /// Whether OpenMCCellAverageProblem should use the displaced mesh
1160 : bool _use_displaced;
1161 :
1162 : /// Mapping from subdomain IDs to which aux variable to read temperature (K) from
1163 : std::map<SubdomainID, std::pair<unsigned int, std::string>> _subdomain_to_temp_vars;
1164 :
1165 : /// Mapping from subdomain IDs to which aux variable to read density (kg/m3) from
1166 : std::map<SubdomainID, std::pair<unsigned int, std::string>> _subdomain_to_density_vars;
1167 :
1168 : /// Mapping from subdomain IDs to the reference density (kg/m3).
1169 : std::map<SubdomainID, Real> _subdomain_to_ref_density;
1170 : };
|
