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 2512 : 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 11197 : virtual const std::map<cellInfo, std::vector<unsigned int>> & cellToElem() const
125 : {
126 11197 : 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 13608 : virtual std::unordered_set<SubdomainID> getCellToElementSub(const cellInfo & info)
135 : {
136 13608 : 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 2297707 : virtual bool hasPointTransformations() const { return _symmetry != nullptr; }
144 :
145 : /**
146 : * Get all the scores added to the tally
147 : * @return scores
148 : */
149 104 : 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 526 : unsigned int getNumScoringTallies(const std::string & score) const
157 : {
158 526 : 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 2248553 : virtual Point transformPoint(const Point & pt) const
243 : {
244 2248553 : 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 25749184 : 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 37346208 : 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 2576128 : 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 5256393 : bool hasDensityFeedback(const cellInfo & cell_info) const
339 : {
340 : std::vector<coupling::CouplingFields> phase = {coupling::density,
341 5256393 : coupling::density_and_temperature};
342 5256393 : return std::find(phase.begin(), phase.end(), cellFeedback(cell_info)) != phase.end();
343 5256393 : }
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 4088024 : bool hasTemperatureFeedback(const cellInfo & cell_info) const
351 : {
352 : std::vector<coupling::CouplingFields> phase = {coupling::temperature,
353 4088024 : coupling::density_and_temperature};
354 4088024 : return std::find(phase.begin(), phase.end(), cellFeedback(cell_info)) != phase.end();
355 4088024 : }
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 1020 : 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 : * Check whether a vector extracted with getParam is empty
429 : * @param[in] vector vector
430 : * @param[in] name name to use for printing error if empty
431 : */
432 : template <typename T>
433 4514 : void checkEmptyVector(const std::vector<T> & vector, const std::string & name) const
434 : {
435 4408 : if (vector.empty())
436 8 : mooseError(name + " cannot be empty!");
437 4506 : }
438 :
439 3834 : int fixedPointIteration() const { return _fixed_point_iteration; }
440 :
441 : /**
442 : * Checks if the problem uses adaptivity or not.
443 : * @return if the problem uses adaptivity.
444 : */
445 2735 : bool hasAdaptivity() const { return _has_adaptivity; }
446 :
447 : /**
448 : * Checks if the problem is using a MoabSkinner or not.
449 : * @return if the problem uses skinning
450 : */
451 133 : bool hasSkinner() const { return _using_skinner; }
452 :
453 : /// Constant flag to indicate that a cell/element was unmapped
454 : static constexpr int32_t UNMAPPED{-1};
455 :
456 : /// Spatial dimension of the Monte Carlo problem
457 : static constexpr int DIMENSION{3};
458 :
459 : /// Get a modifyable non-const reference to the Moose mesh
460 : virtual MooseMesh & getMooseMesh();
461 :
462 : /// Get a modifyable const reference to the Moose mesh
463 : virtual const MooseMesh & getMooseMesh() const;
464 :
465 : /**
466 : * Whether a moving mesh is used
467 : * @return whether the [Mesh] is moving
468 : */
469 : const bool & useDisplaced() const { return _use_displaced; }
470 :
471 : protected:
472 : /**
473 : * A function to re-initialize coupling and apply feedback to the OpenMC problem.
474 : * Applied before OpenMC is executed in either: i) normal Picard iterations, or
475 : * ii) criticality searches. This function performs these operations in the
476 : * following order:
477 : * 1. Updates the OpenMC geometry using the skinner;
478 : * 2. Resets tallies and the cell->element maps to take into account
479 : * mesh/geometry changes;
480 : * 3. Updates the nuclide composition of OpenMC materials;
481 : * 4. Sets cell temperatures and densities;
482 : * 5. Exports OpenMC properties;
483 : * 6. Reinitializes multi-group cross sections to take into account changing
484 : * temperatures and densities.
485 : */
486 : void reinitCouplingAndApplyFeedback();
487 :
488 : /**
489 : * Implement critSearchStep() to re-generate the cell-to-element (and dual)
490 : * mapping. This sends new temperatures and densities to OpenMC from the
491 : * re-mapped elements to ensure the state remains critical under changes to the
492 : * model with feedbacks.
493 : */
494 : virtual void critSearchStep() override;
495 :
496 : /**
497 : * Get the cell level in OpenMC to use for coupling
498 : * @param[in] c point
499 : * @return cell level
500 : */
501 : unsigned int getCellLevel(const Point & c) const;
502 :
503 : /**
504 : * Read the names of the MOOSE variables used for sending feedback into OpenMC
505 : * @param[in] param feedback term to read
506 : * @param[in] default_name default name to use for MOOSE variables holding this field
507 : * @param[out] vars_to_specified_blocks map from MOOSE variable names to the blocks on which they
508 : * are defined
509 : * @param[out] specified_blocks user-specified blocks for feedback
510 : */
511 : void
512 : readBlockVariables(const std::string & param,
513 : const std::string & default_name,
514 : std::map<std::string, std::vector<SubdomainName>> & vars_to_specified_blocks,
515 : std::vector<SubdomainID> & specified_blocks);
516 :
517 : /**
518 : * Whether this cell has an identical fill
519 : * @param[in] cell_info cell index, instance pair
520 : * @return whether this cell has an identical fill
521 : */
522 : bool cellHasIdenticalFill(const cellInfo & cell_info) const;
523 :
524 : /**
525 : * When using the 'identical_cell_fills' feature, this is used to determine the
526 : * contained material cells in each parent cell by applying a uniform shift
527 : * @param[in] cell_info cell index, instance pair
528 : * @return material cells contained within the given cell
529 : */
530 : containedCells shiftCellInstances(const cellInfo & cell_info) const;
531 :
532 : /**
533 : * Whether this cell overlaps with ANY value in the given subdomain set
534 : * @param[in] cell_info cell index, instance pair
535 : * @param[in] id subdomain IDs
536 : * @return whether the cell overlaps with the subdomain
537 : */
538 : bool cellMapsToSubdomain(const cellInfo & cell_info,
539 : const std::unordered_set<SubdomainID> & id) const;
540 :
541 : /**
542 : * Get the first material cell contained in the given cell
543 : * @param[in] cell_info cell index, instance pair
544 : * @return material cell index, instance pair
545 : */
546 : cellInfo firstContainedMaterialCell(const cellInfo & cell_info) const;
547 :
548 : /**
549 : * Get all of the material cells contained within this cell
550 : * @param[in] cell_info cell index, instance pair
551 : * @return all material cells contained in the given cell
552 : */
553 : containedCells containedMaterialCells(const cellInfo & cell_info) const;
554 :
555 : /**
556 : * Delete the OpenMC DAGMC geometry and re-generate the CSG geometry data structures in-place.
557 : */
558 : void updateOpenMCGeometry();
559 :
560 : /**
561 : * Get a list of each material in the problem, sorted by subdomain. This function also checks
562 : * that there is just one OpenMC material in each subdomain, necessary for the DAGMC skinning.
563 : * @return material in each subdomain
564 : */
565 : std::vector<std::string> getMaterialInEachSubdomain() const;
566 :
567 : /**
568 : * Apply transformations and scale point from MOOSE into the OpenMC domain
569 : * @param[in] pt point
570 : * @return transformed point
571 : */
572 : Point transformPointToOpenMC(const Point & pt) const;
573 :
574 : /**
575 : * For geometries with fine-scale details (e.g. TRISO), Cardinal's default settings can
576 : * take a very long time to initialize the problem (but we can't change those defaults
577 : * because they are not 100% applicable all the time). So, we print out a message to
578 : * the user to point them in the right direction if their initialization is taking a
579 : * long time.
580 : * @param[in] start time to use for evaluating whether we've exceeded our limit for printing the
581 : * message
582 : */
583 : void
584 : printTrisoHelp(const std::chrono::time_point<std::chrono::high_resolution_clock> & start) const;
585 :
586 : /**
587 : * Print to the console the names of the auxvariables used for I/O with OpenMC.
588 : * We only print these tables once, upon initialization, because this data does
589 : * not change if re-initializing the spatial mapping for moving-mesh problems,
590 : * adaptive refinement, skinning, etc.
591 : */
592 : void printAuxVariableIO();
593 :
594 : /**
595 : * Get all the material indices within the set of cells
596 : * @param[in] contained_cells set of cells
597 : * @return contained materials
598 : */
599 : std::vector<int32_t> materialsInCells(const containedCells & contained_cells) const;
600 :
601 : /// Loop over the mapped cells, and build a map between subdomains to OpenMC materials
602 : void subdomainsToMaterials();
603 :
604 : /**
605 : * Get a set of all subdomains that have at least 1 element coupled to an OpenMC cell
606 : * @return subdomains with at least 1 element coupled to OpenMC
607 : */
608 : std::set<SubdomainID> coupledSubdomains() const;
609 :
610 : /**
611 : * Gather a vector of values to be summed for each cell
612 : * @param[in] local local values to be summed for the cells
613 : * @param[out] global global mapping of the summed values to the cells
614 : */
615 : template <typename T>
616 : void gatherCellSum(std::vector<T> & local, std::map<cellInfo, T> & global) const;
617 :
618 : /**
619 : * Gather a vector of values to be pushed back to for each cell
620 : * @param[in] local local values to be pushed back for the cells
621 : * @param[in] n_local number of local values contributed to each cell
622 : * @param[out] global global mapping of the pushed back values to the cells
623 : */
624 : template <typename T>
625 : void gatherCellVector(std::vector<T> & local, std::vector<unsigned int> & n_local, std::map<cellInfo, std::vector<T>> & global);
626 :
627 : /**
628 : * Get the feedback which this element provides to OpenMC
629 : * @param[in] elem
630 : * @return coupling phase
631 : */
632 : coupling::CouplingFields elemFeedback(const Elem * elem) const;
633 :
634 : /**
635 : * Read the parameters needed for triggers
636 : * @param[in] params input parameters
637 : */
638 : void getTallyTriggerParameters(const InputParameters & params);
639 :
640 : /**
641 : * Read the block parameters based on user settings
642 : * @param[in] name name of input parameter representing a vector of subdomain names
643 : * @param[in] blocks list of block ids to write
644 : */
645 : void readBlockParameters(const std::string name, std::unordered_set<SubdomainID> & blocks);
646 :
647 : /**
648 : * Cache the material cells contained within each coupling cell;
649 : * depending on user settings, this may attempt to take shortcuts
650 : * by assuming each cell has the same fills
651 : */
652 : void cacheContainedCells();
653 :
654 : /**
655 : * Fill the cached contained cells data structure for a given cell
656 : * @param[in] cell_info cell index, instance pair
657 : * @param[in] hint location hint used to accelerate the search
658 : * @param[out] map contained cell map
659 : */
660 : void setContainedCells(const cellInfo & cell_info,
661 : const Point & hint,
662 : std::map<cellInfo, containedCells> & map);
663 :
664 : /**
665 : * Check that the structure of the contained material cells for two cell matches;
666 : * i.e. this checks that the keys are the same and that the *number* of instances
667 : * of each filling material cell match.
668 : * @param[in] cell_info cell index, instance pair
669 : * @param[in] reference map we want to check against
670 : * @param[in] compare map we want to check
671 : */
672 : void checkContainedCellsStructure(const cellInfo & cell_info,
673 : containedCells & reference,
674 : containedCells & compare) const;
675 :
676 : /**
677 : * Set a minimum order for a volume quadrature rule
678 : * @param[in] volume_order order of the volume quadrature rule
679 : * @param[in] type string type of quadrature rule for printing a console message
680 : */
681 : void setMinimumVolumeQRules(Order & volume_order, const std::string & type);
682 :
683 : /// For keeping the output neat when using verbose
684 : std::string printNewline() const
685 : {
686 : if (_verbose)
687 : return "\n";
688 : else
689 : return "";
690 : }
691 :
692 : /// Loop over the elements in the MOOSE mesh and store the type of feedback applied by each.
693 : void storeElementPhase();
694 :
695 : /**
696 : * Loop over all the OpenMC cells and count the number of MOOSE elements to which the cell
697 : * is mapped based on phase.
698 : */
699 : void getCellMappedPhase();
700 :
701 : /// This function is used to ensure that each OpenMC cell only maps to a single phase
702 : void checkCellMappedPhase();
703 :
704 : /// Loop over all the OpenMC cells and get the element subdomain IDs that map to each cell
705 : void getCellMappedSubdomains();
706 :
707 : /**
708 : * Loop over all the OpenMC cells and compute the volume of the MOOSE elements that each
709 : * cell maps to
710 : */
711 : void computeCellMappedVolumes();
712 :
713 : /// Set up the mapping from MOOSE elements to OpenMC cells
714 : void initializeElementToCellMapping();
715 :
716 : /// Populate maps of MOOSE elements to OpenMC cells
717 : void mapElemsToCells();
718 :
719 : /**
720 : * A function which validates local tallies. This is done to ensure that at least one of the
721 : * tallies contains a heating score when running in eigenvalue mode. This must be done outside
722 : * of the constructor as tallies are added from an external system.
723 : */
724 : void validateLocalTallies();
725 :
726 : /// Add OpenMC tallies to facilitate the coupling
727 : void initializeTallies();
728 :
729 : /**
730 : * Reset any tallies previously added by Cardinal, by deleting them from OpenMC.
731 : * Also delete any mesh filters and meshes added to OpenMC for mesh filters.
732 : */
733 : void resetTallies();
734 :
735 : /// Find the material filling each cell which receives density feedback
736 : void getMaterialFills();
737 :
738 : /**
739 : * Get one point inside each cell, for accelerating the particle search routine.
740 : * This function will get the centroid of the first global element in the lowest
741 : * rank in the cell.
742 : */
743 : void getPointInCell();
744 :
745 : /**
746 : * Compute the product of volume with a field across ranks and sum into a global map
747 : * @param[in] var_num variable to weight with volume, mapped by subdomain ID
748 : * @param[in] phase phases to compute the operation for
749 : * @return volume-weighted field for each cell, in a global sense
750 : */
751 : std::map<cellInfo, Real> computeVolumeWeightedCellInput(
752 : const std::map<SubdomainID, std::pair<unsigned int, std::string>> & var_num,
753 : const std::vector<coupling::CouplingFields> * phase) const;
754 :
755 : /**
756 : * Send temperature from MOOSE to OpenMC by computing a volume average
757 : * and applying a single temperature per OpenMC cell
758 : */
759 : void sendTemperatureToOpenMC() const;
760 :
761 : /**
762 : * Send density from MOOSE to OpenMC by computing a volume average
763 : * and applying a single density per OpenMC cell.
764 : */
765 : void sendDensityToOpenMC() const;
766 :
767 : /**
768 : * Check if a mapped location is in the outer universe of a lattice
769 : * @param[in] level lattice level
770 : * @return whether the location is in the outer universe
771 : */
772 : void latticeOuterCheck(const Point & c, int level) const;
773 :
774 : /**
775 : * Report an error for a mapped location in an outer universe of a lattice
776 : * @param[in] c Mapped location
777 : * @param[in] level level of the mapped cell
778 : */
779 : void latticeOuterError(const Point & c, int level) const;
780 :
781 : /**
782 : * Find the OpenMC cell at a given point in space
783 : * @param[in] point point
784 : * @return whether OpenMC reported an error
785 : */
786 : bool findCell(const Point & point);
787 :
788 : /**
789 : * Checks that the contained material cells exactly match between a reference obtained
790 : * by calling openmc::Cell::get_contained_cells for each cell and a shortcut
791 : * approach that assumes all identical cells (which aren't simply just material fills)
792 : * has exactly the same contained material cells.
793 : * @param[in] reference reference map to compare against
794 : * @param[in] compare shortcut map to compare
795 : */
796 : void compareContainedCells(std::map<cellInfo, containedCells> & reference,
797 : std::map<cellInfo, containedCells> & compare) const;
798 :
799 : NumericVector<Number> & _serialized_solution;
800 :
801 : /**
802 : * Return all IDs of all Cardinal-mapped Tallies
803 : * @return all Cardinal-mapped Tally IDs
804 : */
805 : virtual std::vector<int32_t> getMappedTallyIDs() const override;
806 :
807 : /**
808 : * Whether to automatically compute the mapping of OpenMC cell IDs and
809 : * instances to the [Mesh].
810 : */
811 : const bool & _output_cell_mapping;
812 :
813 : /**
814 : * Where to get the initial OpenMC temperatures and densities from;
815 : * can be either hdf5 (from a properties.h5 file), xml (whatever is already
816 : * set in the XML files), or moose (meaning whatever ICs are set on the 'temperature_variables'
817 : * and 'density_variables'
818 : */
819 : const coupling::OpenMCInitialCondition _initial_condition;
820 :
821 : /// Type of relaxation to apply to the OpenMC tallies
822 : const relaxation::RelaxationEnum _relaxation;
823 :
824 : /**
825 : * Type of trigger to apply to k eigenvalue to indicate when
826 : * the simulation is complete. These can be used to on-the-fly adjust the number
827 : * of active batches in order to reach some desired criteria (which is specified
828 : * by this parameter).
829 : */
830 : const trigger::TallyTriggerTypeEnum _k_trigger;
831 :
832 : /**
833 : * Coordinate level in the OpenMC domain to use for mapping cells to mesh.
834 : * When using 'lowest_cell_level', this parameter indicates that the lowest
835 : * cell level is used, up until _cell_level.
836 : */
837 : unsigned int _cell_level;
838 :
839 : /**
840 : * Whether OpenMC properties (temperature and density) should be exported
841 : * after being updated in syncSolutions.
842 : */
843 : const bool & _export_properties;
844 :
845 : /// Whether or not the problem uses a skinner to regenerate the OpenMC geometry.
846 : const bool _using_skinner;
847 :
848 : /**
849 : * When the mesh changes during the simulation (either from adaptive mesh refinement
850 : * or deformation), the mapping from OpenMC cells to the [Mesh] must be re-established
851 : * after each OpenMC run.
852 : */
853 : bool _need_to_reinit_coupling;
854 :
855 : /**
856 : * If known a priori by the user, whether the tally cells (which are not simply material
857 : * fills) have EXACTLY the same contained material cells. This is a big optimization for
858 : * TRISO problems in setting up homogenized temperature/density feedback to OpenMC.
859 : *
860 : * The concept can best be explained with a pebble bed reactor.
861 : * If every pebble is filled with an identical TRISO universe, then the material fills
862 : * in each pebble are identical to one another except for a constant offset. This idea
863 : * can be used to then skip all but the first two openmc::Cell::get_contained_cells
864 : * calls (which are required in order to figure out the pattern by which pebble N is
865 : * incremented relative to pebble 1).
866 : *
867 : * When using this parameter, we HIGHLY recommend setting 'check_identical_cell_fills =
868 : * true' the first time you run your model. This will figure out the material cell fills using a
869 : * method that calls openmc::Cell::get_contained_cells for every tally cell, i.e. without assuming
870 : * anything about repeated structure in your OpenMC model. Setting 'identical_cell_fills'
871 : * without also setting 'check_identical_cell_fills = true' may result in SILENT
872 : * errors!!! So it is essential to be sure you've removed any error sources before you turn the
873 : * error check off to actually leverage the speedup.
874 : *
875 : * Note: for any tally cells that are just filled with a material, we use the approach
876 : * where openmc::Cell::get_contained_cells is called in full.
877 : *
878 : * This optimization will not work (and 'check_identical_cell_fills = true' *will*
879 : * catch these) for:
880 : * - any situation where tallied, non-material-fill pebbles have different fills
881 : * (such as if you have different TRISO lattices in each pebble)
882 : * - any situation where there is a "gap" in the incrementing of the material fill
883 : * instances (such as if pebble 89 does not map to 'tally_blocks', then the instance
884 : * shift for pebble 90 relative to pebble 1 is 89, when it should have been 90).
885 : */
886 : const bool _has_identical_cell_fills;
887 :
888 : /**
889 : * Whether we should rigorously check that each tally cell has identical fills;
890 : * this is SLOW for large TRISO problems, but is essential to ensure the accuracy of
891 : * 'identical_cell_fills'. Please set 'check_identical_cell_fills' to 'true' at least
892 : * once before running production cases to be sure the optimization can be applied.
893 : */
894 : const bool & _check_identical_cell_fills;
895 :
896 : /**
897 : * Whether it can be assumed that all of the tallies (both those set by the user
898 : * in the XML file, as well as those created automatically by Cardinal) are
899 : * spatially separate. This means that once a particle scores to one tally bin, it wouldn't
900 : * score to ANY other tally bins. This can dramatically increase tracking rates
901 : * for problems with many tallies.
902 : */
903 : const bool & _assume_separate_tallies;
904 :
905 : /**
906 : * Whether to map density according to each individual OpenMC cell (in which case an
907 : * error is thrown if you don't have a unique material in each cell) or by material.
908 : */
909 : bool _map_density_by_cell;
910 :
911 : /**
912 : * Whether the problem has density feedback blocks specified; note that this is NOT necessarily
913 : * indicative that the mapping was successful in finding any cells corresponding to those blocks
914 : */
915 : const bool _specified_density_feedback;
916 :
917 : /**
918 : * Whether the problem has temperature feedback blocks specified; note that this is NOT
919 : * necessarily indicative that the mapping was successful in finding any cells corresponding to
920 : * those blocks
921 : */
922 : const bool _specified_temperature_feedback;
923 :
924 : /// Whether any cell tallies exist.
925 : bool _has_cell_tallies = false;
926 :
927 : /// Whether any spatial mapping from OpenMC's cells to the mesh is needed
928 : bool _needs_to_map_cells;
929 :
930 : /**
931 : * A map of the filter objects created by the [Problem/Filters] block. The key for each filter is
932 : * it's corresponding MOOSE name to allow tallies to look up filters.
933 : */
934 : std::map<std::string, std::shared_ptr<FilterBase>> _filters;
935 :
936 : /// A vector of the tally objects created by the [Problem/Tallies] block.
937 : std::vector<std::shared_ptr<TallyBase>> _local_tallies;
938 :
939 : /// A list of all of the scores contained by the local tallies added in the [Tallies] block.
940 : std::vector<std::string> _all_tally_scores;
941 :
942 : /// Number of tallies scoring a particular score.
943 : std::map<std::string, unsigned int> _score_count;
944 :
945 : /// A vector of auxvariable ids added by the [Tallies] block.
946 : std::vector<std::vector<unsigned int>> _tally_var_ids;
947 :
948 : /// A vector of external (output-based) auxvariable ids added by the [Tallies] block.
949 : std::vector<std::vector<std::vector<unsigned int>>> _tally_ext_var_ids;
950 :
951 : /// Blocks in MOOSE mesh that provide density feedback
952 : std::vector<SubdomainID> _density_blocks;
953 :
954 : /// Blocks in MOOSE mesh that provide temperature feedback
955 : std::vector<SubdomainID> _temp_blocks;
956 :
957 : /// Blocks for which the cell fills are identical
958 : std::unordered_set<SubdomainID> _identical_cell_fill_blocks;
959 :
960 : /// Mapping of MOOSE elements to the OpenMC cell they map to (if any)
961 : std::vector<cellInfo> _elem_to_cell{};
962 :
963 : /// Phase of each cell
964 : std::map<cellInfo, coupling::CouplingFields> _cell_phase;
965 :
966 : /// Number of elements in the MOOSE mesh that exclusively provide density feedback
967 : int _n_moose_density_elems;
968 :
969 : /// Number of elements in the MOOSE mesh that exclusively provide temperature feedback
970 : int _n_moose_temp_elems;
971 :
972 : /// Number of elements in the MOOSE mesh which provide temperature+density feedback
973 : int _n_moose_temp_density_elems;
974 :
975 : /// Number of no-coupling elements in the MOOSE mesh
976 : int _n_moose_none_elems;
977 :
978 : /**
979 : * Number of MOOSE elements that exclusively provide temperature feedback,
980 : * and which successfully mapped to OpenMC cells
981 : */
982 : int _n_mapped_temp_elems;
983 :
984 : /**
985 : * Number of MOOSE elements that exclusively provide density feedback,
986 : * and which successfully mapped to OpenMC cells
987 : */
988 : int _n_mapped_density_elems;
989 :
990 : /**
991 : * Number of MOOSE elements that provide temperature+density feedback,
992 : * and which successfully mapped to OpenMC cells
993 : */
994 : int _n_mapped_temp_density_elems;
995 :
996 : /// Number of no-coupling elements mapped to OpenMC cells
997 : int _n_mapped_none_elems;
998 :
999 : /// Total volume of uncoupled MOOSE mesh elements
1000 : Real _uncoupled_volume;
1001 :
1002 : /// Whether non-material cells are mapped
1003 : bool _material_cells_only{true};
1004 :
1005 : /// Mapping of OpenMC cell indices to a vector of MOOSE element IDs
1006 : std::map<cellInfo, std::vector<unsigned int>> _cell_to_elem;
1007 :
1008 : /// Mapping of OpenMC cell indices to a vector of MOOSE element IDs, on each local rank
1009 : std::map<cellInfo, std::vector<unsigned int>> _local_cell_to_elem;
1010 :
1011 : /// Mapping of OpenMC cell indices to the subdomain IDs each maps to
1012 : std::map<cellInfo, std::unordered_set<SubdomainID>> _cell_to_elem_subdomain;
1013 :
1014 : /// Mapping of elem subdomains to materials
1015 : std::map<SubdomainID, std::set<int32_t>> _subdomain_to_material;
1016 :
1017 : /**
1018 : * A point inside the cell, taken simply as the centroid of the first global
1019 : * element inside the cell. This is stored to accelerate the particle search.
1020 : */
1021 : std::map<cellInfo, Point> _cell_to_point;
1022 :
1023 : /**
1024 : * Volume associated with the mapped element space for each OpenMC cell; the unit
1025 : * for this volume is whatever is used in the [Mesh] block
1026 : */
1027 : std::map<cellInfo, Real> _cell_to_elem_volume;
1028 :
1029 : /**
1030 : * Volume associated with the actual OpenMC cell, computed by an optional
1031 : * OpenMCVolumeCalculation user object
1032 : */
1033 : std::map<cellInfo, Real> _cell_volume;
1034 :
1035 : /**
1036 : * Material filling each cell to receive density feedback. We enforce that these
1037 : * cells are filled with a material (cannot be filled with a lattice or universe).
1038 : */
1039 : std::map<cellInfo, int32_t> _cell_to_material;
1040 :
1041 : /**
1042 : * Material-type cells contained within a cell; this is only populated if a cell
1043 : * is NOT indicated as having an identical fill
1044 : */
1045 : std::map<cellInfo, containedCells> _cell_to_contained_material_cells;
1046 :
1047 : /// Number of material-type cells contained within a cell
1048 : std::map<cellInfo, int32_t> _cell_to_n_contained;
1049 :
1050 : /// Whether the present transfer is the first transfer
1051 : static bool _first_transfer;
1052 :
1053 : /// Whether the diagnostic tables on initialization have already been printed
1054 : static bool _printed_initial;
1055 :
1056 : /// Whether a warning has already been printed about very long setup times (for TRISOs)
1057 : static bool _printed_triso_warning;
1058 :
1059 : /// Dummy particle to reduce number of allocations of particles for cell lookup routines
1060 : openmc::Particle _particle;
1061 :
1062 : /// Number of particles simulated in the first iteration in Dufek-Gudowski relaxation
1063 : unsigned int _n_particles_1;
1064 :
1065 : /// Mapping from temperature variable name to the subdomains on which to read it from
1066 : std::map<std::string, std::vector<SubdomainName>> _temp_vars_to_blocks;
1067 :
1068 : /// Mapping from density variable name to the subdomains on which to read it from
1069 : std::map<std::string, std::vector<SubdomainName>> _density_vars_to_blocks;
1070 :
1071 : /// Optional volume calculation for cells which map to MOOSE
1072 : OpenMCVolumeCalculation * _volume_calc;
1073 :
1074 : /// Userobject that maps from a partial-symmetry OpenMC model to a whole-domain [Mesh]
1075 : const SymmetryPointGenerator * _symmetry;
1076 :
1077 : /// Number of temperature-only feedback elements in each mapped OpenMC cell (global)
1078 : std::map<cellInfo, int> _n_temp;
1079 :
1080 : /// Number of density-only feedback elements in each mapped OpenMC cell (global)
1081 : std::map<cellInfo, int> _n_rho;
1082 :
1083 : /// Number of temperature+density feedback elements in each mapped OpenMC cell (global)
1084 : std::map<cellInfo, int> _n_temp_rho;
1085 :
1086 : /// Number of none elements in each mapped OpenMC cell (global)
1087 : std::map<cellInfo, int> _n_none;
1088 :
1089 : /// The tally to be used for normalizing all other tallies when running an eigenvalue calculation.
1090 : std::shared_ptr<TallyBase> _source_rate_norm_tally;
1091 :
1092 : /// The score index into "_source_rate_norm_tally".
1093 : int _source_rate_score = -1;
1094 :
1095 : #ifdef ENABLE_DAGMC
1096 : /// Optional skinner to re-generate the OpenMC geometry on-the-fly for DAGMC models
1097 : MoabSkinner * _skinner = nullptr;
1098 :
1099 : /// Pointer to DAGMC
1100 : std::shared_ptr<moab::DagMC> _dagmc = nullptr;
1101 : #endif
1102 :
1103 : /// Total number of unique OpenMC cell IDs + instances combinations
1104 : long unsigned int _n_openmc_cells;
1105 :
1106 : /// ID of the OpenMC universe corresponding to the DAGMC universe
1107 : int32_t _dagmc_universe_id;
1108 :
1109 : /// Whether the DAGMC universe is the root universe or not.
1110 : bool _dagmc_root_universe = true;
1111 :
1112 : /// ID of the OpenMC cell corresponding to the cell which uses the DAGMC universe as a fill.
1113 : int32_t _cell_using_dagmc_universe_id;
1114 :
1115 : /// The number of OpenMC surfaces before skinning occurs. This is required to properly reinitialize
1116 : /// the CSG geometry contained in the OpenMC model.
1117 : const int32_t _initial_num_openmc_surfaces;
1118 :
1119 : /// Conversion rate from eV to Joule
1120 : static constexpr Real EV_TO_JOULE = 1.6022e-19;
1121 :
1122 : /// Tolerance for setting zero tally
1123 : static constexpr Real ZERO_TALLY_THRESHOLD = 1e-12;
1124 :
1125 : private:
1126 : /**
1127 : * Update the number of particles according to the Dufek-Gudowski relaxation scheme
1128 : */
1129 : void dufekGudowskiParticleUpdate();
1130 :
1131 : /// Flattened cell IDs collected after parallel communication
1132 : std::vector<int32_t> _flattened_ids;
1133 :
1134 : /// Flattened cell instancess collected after parallel communication
1135 : std::vector<int32_t> _flattened_instances;
1136 :
1137 : /// Offsets for each cell instance in an identically-repeated universe
1138 : containedCells _instance_offsets;
1139 :
1140 : /// Offset for each cell relative to the first identical-fill cell
1141 : std::map<cellInfo, int32_t> _n_offset;
1142 :
1143 : /// First identical-fill cell
1144 : cellInfo _first_identical_cell;
1145 :
1146 : /// Materials in the first identical-fill cell
1147 : std::vector<int32_t> _first_identical_cell_materials;
1148 :
1149 : /// Whether OpenMCCellAverageProblem should use the displaced mesh
1150 : bool _use_displaced;
1151 :
1152 : /// Mapping from subdomain IDs to which aux variable to read temperature (K) from
1153 : std::map<SubdomainID, std::pair<unsigned int, std::string>> _subdomain_to_temp_vars;
1154 :
1155 : /// Mapping from subdomain IDs to which aux variable to read density (kg/m3) from
1156 : std::map<SubdomainID, std::pair<unsigned int, std::string>> _subdomain_to_density_vars;
1157 : };
|
