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 : #ifdef ENABLE_OPENMC_COUPLING
20 : #include "TallyBase.h"
21 :
22 : #include "OpenMCCellAverageProblem.h"
23 : #include "UserErrorChecking.h"
24 : #include "AuxiliarySystem.h"
25 : #include "FilterBase.h"
26 :
27 : #include "AngularLegendreFilter.h"
28 : #include "EnergyOutFilter.h"
29 : #include "DelayedGroupFilter.h"
30 :
31 : #include "openmc/tallies/filter_universe.h"
32 : #include "openmc/geometry.h"
33 : #include "openmc/settings.h"
34 : #include "openmc/universe.h"
35 :
36 : InputParameters
37 5378 : TallyBase::validParams()
38 : {
39 5378 : auto params = MooseObject::validParams();
40 10756 : params.addParam<MultiMooseEnum>(
41 : "score",
42 10756 : getTallyScoreEnum(),
43 : "Score(s) to use in the OpenMC tallies. If not specified, defaults to 'kappa_fission'");
44 10756 : params.addParam<MooseEnum>(
45 10756 : "estimator", getTallyEstimatorEnum(), "Type of tally estimator to use in OpenMC");
46 10756 : params.addParam<std::vector<std::string>>(
47 : "name",
48 : "Auxiliary variable name(s) to use for OpenMC tallies. "
49 : "If not specified, defaults to the names of the scores");
50 :
51 10756 : params.addParam<std::vector<SubdomainName>>(
52 : "block",
53 : "Subdomains for which to add tallies in OpenMC. If not provided, "
54 : "tallies will be applied over the entire domain corresponding to the [Mesh] block.");
55 10756 : params.addParam<std::vector<SubdomainName>>("blocks",
56 : "This parameter is deprecated, use 'block' instead!");
57 :
58 5378 : MultiMooseEnum tally_trigger("rel_err none");
59 10756 : params.addParam<MultiMooseEnum>(
60 : "trigger",
61 : tally_trigger,
62 : "Trigger criterion to determine when OpenMC simulation is complete "
63 : "based on tallies. If multiple scores are specified in 'score, "
64 : "this same trigger is applied to all scores.");
65 10756 : params.addRangeCheckedParam<std::vector<Real>>(
66 : "trigger_threshold", "trigger_threshold > 0", "Threshold for the tally trigger");
67 10756 : params.addParam<std::vector<bool>>(
68 : "trigger_ignore_zeros",
69 : {false},
70 : "Whether tally bins with zero scores are ignored when computing the tally trigger. If only "
71 : "one "
72 : "value of 'trigger_ignore_zeros' is provided, that value is applied to all tally scores.");
73 :
74 : MultiMooseEnum openmc_outputs(
75 5378 : "unrelaxed_tally_std_dev unrelaxed_tally_rel_error unrelaxed_tally");
76 10756 : params.addParam<MultiMooseEnum>(
77 : "output",
78 : openmc_outputs,
79 : "UNRELAXED field(s) to output from OpenMC for each tally score. "
80 : "unrelaxed_tally_std_dev will write the standard deviation of "
81 : "each tally into auxiliary variables "
82 : "named *_std_dev. unrelaxed_tally_rel_error will write the "
83 : "relative standard deviation (unrelaxed_tally_std_dev / unrelaxed_tally) "
84 : "of each tally into auxiliary variables named *_rel_error. "
85 : "unrelaxed_tally will write the raw unrelaxed tally into auxiliary "
86 : "variables named *_raw (replace * with 'name').");
87 :
88 10756 : params.addParam<std::vector<std::string>>("filters", "External filters to add to this tally.");
89 :
90 10756 : params.addParam<bool>("check_tally_sum",
91 : "Whether to check consistency between the local tallies "
92 : "with a global tally sum. This will require that the "
93 : "integral of the local tally matches a tally with no filters "
94 : "(defined over the entire phase space).");
95 16134 : params.addRangeCheckedParam<Real>("tally_sum_tol",
96 10756 : 1e-6,
97 : "tally_sum_tol > 0",
98 : "Relative tolerance above which to throw an error message that "
99 : "the local tally does not preserve the global tally");
100 10756 : params.addParam<bool>(
101 : "normalize_by_global_tally",
102 10756 : true,
103 : "Whether to normalize local tallies by a global tally (true) or else by the sum "
104 : "of the local tally (false)");
105 :
106 5378 : params.addPrivateParam<OpenMCCellAverageProblem *>("_openmc_problem");
107 :
108 5378 : params.registerBase("Tally");
109 5378 : params.registerSystemAttributeName("Tally");
110 :
111 5378 : return params;
112 5378 : }
113 :
114 2839 : TallyBase::TallyBase(const InputParameters & parameters)
115 : : MooseObject(parameters),
116 2839 : _openmc_problem(*getParam<OpenMCCellAverageProblem *>("_openmc_problem")),
117 2839 : _mesh(_openmc_problem.mesh()),
118 2839 : _aux(_openmc_problem.getAuxiliarySystem()),
119 3067 : _tally_trigger(isParamValid("trigger") ? &getParam<MultiMooseEnum>("trigger") : nullptr),
120 8517 : _trigger_ignore_zeros(getParam<std::vector<bool>>("trigger_ignore_zeros")),
121 2839 : _normalize_by_global(_openmc_problem.runMode() == openmc::RunMode::FIXED_SOURCE
122 2839 : ? false
123 5532 : : getParam<bool>("normalize_by_global_tally")),
124 2839 : _check_tally_sum(isParamValid("check_tally_sum")
125 6349 : ? getParam<bool>("check_tally_sum")
126 2168 : : (_openmc_problem.runMode() == openmc::RunMode::FIXED_SOURCE
127 2168 : ? true
128 2074 : : _normalize_by_global)),
129 5678 : _tally_sum_tol(getParam<Real>("tally_sum_tol")),
130 2839 : _needs_global_tally(_check_tally_sum || _normalize_by_global),
131 5678 : _renames_tally_vars(isParamValid("name")),
132 5678 : _has_outputs(isParamValid("output")),
133 8517 : _is_adaptive(_openmc_problem.hasAdaptivity())
134 : {
135 5678 : if (isParamValid("score"))
136 : {
137 3428 : const auto & scores = getParam<MultiMooseEnum>("score");
138 3948 : for (const auto & score : scores)
139 4468 : _tally_score.push_back(_openmc_problem.enumToTallyScore(score));
140 : }
141 : else
142 2250 : _tally_score = {"kappa-fission"};
143 :
144 2839 : if (_openmc_problem.runRandomRay())
145 : {
146 : const auto has_flux =
147 102 : std::find(_tally_score.begin(), _tally_score.end(), "flux") != _tally_score.end();
148 102 : if (openmc::FlatSourceDomain::volume_normalized_flux_tallies_ && has_flux)
149 2 : mooseError("Cardinal volume-normalizes flux tallies with volumes computed from the "
150 : "mesh elements, and so normalizing flux tallies with the random ray "
151 : "volume estimator will result in dividing by volume twice. Please set "
152 : "openmc.Settings.random_ray['volume_estimator'] = False in your OpenMC "
153 : "model.");
154 :
155 : bool found_invalid_rr_score = false;
156 : std::string invalid_scores;
157 302 : for (const auto & score : _tally_score)
158 : {
159 202 : if (!_openmc_problem.validRandomRayScore(score))
160 : {
161 52 : invalid_scores += score + ", ";
162 : found_invalid_rr_score = true;
163 : }
164 : }
165 100 : if (found_invalid_rr_score)
166 : {
167 2 : invalid_scores.erase(invalid_scores.length() - 2, 2);
168 2 : paramError("score",
169 2 : "OpenMC's random ray solver currently doesn't support the following scores: " +
170 0 : invalid_scores + ". Please remove these scores.");
171 : }
172 : }
173 :
174 : const bool heating =
175 2835 : std::find(_tally_score.begin(), _tally_score.end(), "heating") != _tally_score.end();
176 : const bool nu_scatter =
177 2835 : std::find(_tally_score.begin(), _tally_score.end(), "nu-scatter") != _tally_score.end();
178 :
179 : // Random ray mode requires tracklength estimators.
180 2835 : if (_openmc_problem.runRandomRay())
181 : {
182 : // Error if the user tries to set something other than tracklength for random ray mode
183 196 : if (isParamValid("estimator"))
184 : {
185 2 : auto estimator = getParam<MooseEnum>("estimator").getEnum<tally::TallyEstimatorEnum>();
186 2 : if (estimator != tally::tracklength)
187 2 : paramError("estimator",
188 : "The random ray solver in OpenMC requires that tallies use "
189 : "tracklength estimators! Please set 'estimator' to 'tracklength'.");
190 : }
191 96 : _estimator = openmc::TallyEstimator::TRACKLENGTH;
192 : }
193 : else
194 : {
195 5474 : if (isParamValid("estimator"))
196 : {
197 310 : auto estimator = getParam<MooseEnum>("estimator").getEnum<tally::TallyEstimatorEnum>();
198 : // Photon heating tallies cannot use tracklength estimators.
199 310 : if (estimator == tally::tracklength && openmc::settings::photon_transport && heating)
200 2 : paramError("estimator",
201 : "Tracklength estimators are currently incompatible with photon transport and "
202 : "heating scores! For more information: https://tinyurl.com/3wre3kwt");
203 :
204 308 : if (estimator != tally::analog && nu_scatter)
205 2 : paramError("estimator", "Non-analog estimators are not supported for nu_scatter scores!");
206 :
207 306 : _estimator = _openmc_problem.tallyEstimator(estimator);
208 : }
209 : else
210 : {
211 : // Set a default of tracklength. This must be overridden in derived classes that use different
212 : // spatial filters (e.g. unstructured mesh tallies).
213 2427 : _estimator = openmc::TallyEstimator::TRACKLENGTH;
214 :
215 : // Heating tallies in photon transport and nu_scatter tallies require collision/analog scores.
216 2427 : if (nu_scatter && !(heating && openmc::settings::photon_transport))
217 8 : _estimator = openmc::TallyEstimator::ANALOG;
218 2419 : else if (nu_scatter && heating && openmc::settings::photon_transport)
219 2 : paramError(
220 : "estimator",
221 : "A single tally cannot score both nu_scatter and heating when photon transport is "
222 : "enabled, as both scores require different estimators. Consider adding one tally "
223 : "which scores nu_scatter (with an analog estimator), and a second tally that scores "
224 : "heating (with a collision estimator).");
225 :
226 2425 : if (heating && openmc::settings::photon_transport)
227 8 : _estimator = openmc::TallyEstimator::COLLISION;
228 : }
229 : }
230 :
231 2827 : if (heating && !openmc::settings::photon_transport)
232 170 : mooseWarning(
233 : "When using the 'heating' score with photon transport disabled, energy deposition\n"
234 : "from photons is neglected unless you specifically ran NJOY to produce MT=301 with\n"
235 : "photon energy deposited locally (not true for any pre-packaged OpenMC data libraries\n"
236 : "on openmc.org).\n\n"
237 : "If you did NOT specifically run NJOY yourself with this customization, we recommend\n"
238 : "using the 'heating_local' score instead, which will capture photon energy deposition.\n"
239 : "Otherwise, you will underpredict the true energy deposition.");
240 :
241 8481 : if (isParamValid("trigger") != isParamValid("trigger_threshold"))
242 2 : paramError("trigger",
243 : "You must either specify none or both of 'trigger' and "
244 : "'trigger_threshold'. You have specified only one.");
245 :
246 2825 : if (_tally_trigger)
247 : {
248 224 : checkRequiredParam(parameters, "trigger_threshold", "using tally triggers");
249 224 : _tally_trigger_threshold = getParam<std::vector<Real>>("trigger_threshold");
250 :
251 112 : if (_tally_trigger->size() != _tally_score.size())
252 2 : paramError("trigger",
253 2 : "'trigger' (size " + std::to_string(_tally_trigger->size()) +
254 2 : ") must have the same length as 'score' (size " +
255 2 : std::to_string(_tally_score.size()) + ")");
256 :
257 110 : if (_tally_trigger_threshold.size() != _tally_score.size())
258 2 : paramError("trigger_threshold",
259 2 : "'trigger_threshold' (size " + std::to_string(_tally_trigger_threshold.size()) +
260 2 : ") must have the same length as 'score' (size " +
261 2 : std::to_string(_tally_score.size()) + ")");
262 :
263 108 : if (_trigger_ignore_zeros.size() > 1)
264 : {
265 2 : if (_tally_score.size() != _trigger_ignore_zeros.size())
266 2 : paramError("trigger_ignore_zeros",
267 2 : "'trigger_ignore_zeros' (size " + std::to_string(_trigger_ignore_zeros.size()) +
268 2 : ") must have the same length as 'score' (size " +
269 2 : std::to_string(_tally_score.size()) + ")");
270 : }
271 106 : else if (_trigger_ignore_zeros.size() == 1)
272 104 : _trigger_ignore_zeros.resize(_tally_score.size(), _trigger_ignore_zeros[0]);
273 :
274 210 : _openmc_problem.checkEmptyVector(_trigger_ignore_zeros, "trigger_ignore_zeros");
275 : }
276 :
277 : // Fetch the filters required by this tally. Error if the filter hasn't been added yet.
278 5634 : if (isParamValid("filters"))
279 : {
280 2394 : for (const auto & filter_name : getParam<std::vector<std::string>>("filters"))
281 : {
282 2110 : if (!_openmc_problem.hasFilter(filter_name))
283 4 : paramError("filters", "Filter with the name " + filter_name + " does not exist!");
284 :
285 1054 : _ext_filters.push_back(_openmc_problem.getFilter(filter_name));
286 : }
287 : }
288 :
289 : // Check the estimator to make sure it doesn't conflict with certain filters.
290 3863 : for (auto & f : _ext_filters)
291 : {
292 1054 : if ((dynamic_cast<AngularLegendreFilter *>(f.get()) ||
293 1054 : dynamic_cast<EnergyOutFilter *>(f.get())) &&
294 186 : _estimator != openmc::TallyEstimator::ANALOG)
295 8 : paramError("estimator",
296 4 : "The filter " + f->name() +
297 : " requires an analog estimator! Please ensure 'estimator' is set to analog.");
298 :
299 1050 : if (dynamic_cast<DelayedGroupFilter *>(f.get()))
300 26 : for (const auto & s : _tally_score)
301 18 : if (s != "delayed-nu-fission" && s != "decay-rate")
302 4 : paramError("score",
303 2 : "The filter " + f->name() +
304 : " can only be used with delayed_nu_fission and decay_rate scores!");
305 : }
306 :
307 5618 : if (isParamValid("name"))
308 1755 : _tally_name = getParam<std::vector<std::string>>("name");
309 : else
310 : {
311 4798 : for (auto score : _tally_score)
312 : {
313 : std::replace(score.begin(), score.end(), '-', '_');
314 2574 : _tally_name.push_back(score);
315 : }
316 : }
317 :
318 2809 : if (_has_outputs)
319 : {
320 : // names of output are appended to ends of 'name'
321 859 : for (const auto & o : getParam<MultiMooseEnum>("output"))
322 : {
323 : std::string name = o;
324 :
325 337 : if (o == "UNRELAXED_TALLY_STD_DEV")
326 320 : _output_name.push_back("std_dev");
327 177 : else if (o == "UNRELAXED_TALLY_REL_ERROR")
328 272 : _output_name.push_back("rel_error");
329 41 : else if (o == "UNRELAXED_TALLY")
330 82 : _output_name.push_back("raw");
331 : else
332 0 : mooseError("Unhandled OutputEnum in OpenMCCellAverageProblem!");
333 : }
334 : }
335 :
336 2809 : if (_tally_name.size() != _tally_score.size())
337 2 : paramError("name", "'name' must be the same length as 'score'!");
338 :
339 : // Modify the variable names so they take into account the bins in the external filters.
340 2807 : auto all_var_names = _tally_name;
341 3855 : for (const auto & filter : _ext_filters)
342 : {
343 : std::vector<std::string> n;
344 2900 : for (unsigned int i = 0; i < all_var_names.size(); ++i)
345 5006 : for (unsigned int j = 0; j < filter->numBins(); ++j)
346 6308 : n.push_back(all_var_names[i] + "_" + filter->binName(j));
347 :
348 1048 : all_var_names = n;
349 :
350 1048 : _num_ext_filter_bins *= filter->numBins();
351 1048 : }
352 2807 : _tally_name = all_var_names;
353 :
354 : // A map of external filter bins to skip when computing sums and means for normalization.
355 2807 : std::vector<bool> skip{false};
356 3855 : for (const auto & filter : _ext_filters)
357 : {
358 : std::vector<bool> s;
359 2678 : for (unsigned int i = 0; i < skip.size(); ++i)
360 4424 : for (unsigned int j = 0; j < filter->numBins(); ++j)
361 3218 : s.push_back(skip[i] || filter->skipBin(j));
362 :
363 1048 : skip = s;
364 : }
365 2807 : _ext_bins_to_skip = skip;
366 :
367 5614 : if (isParamSetByUser("blocks"))
368 0 : mooseError("This parameter is deprecated, use 'block' instead!");
369 :
370 5614 : if (isParamValid("block"))
371 : {
372 5070 : auto block_names = getParam<std::vector<SubdomainName>>("block");
373 1690 : if (block_names.empty())
374 2 : paramError("block", "Subdomain names must be provided if using 'block'!");
375 :
376 1688 : auto block_ids = _openmc_problem.getMooseMesh().getSubdomainIDs(block_names);
377 1688 : std::copy(
378 1688 : block_ids.begin(), block_ids.end(), std::inserter(_tally_blocks, _tally_blocks.end()));
379 :
380 : // Check to make sure all of the blocks are in the mesh.
381 1688 : const auto & subdomains = _openmc_problem.getMooseMesh().meshSubdomains();
382 4349 : for (std::size_t b = 0; b < block_names.size(); ++b)
383 2663 : if (subdomains.find(block_ids[b]) == subdomains.end())
384 4 : paramError("block",
385 2 : "Block '" + block_names[b] + "' specified in 'block' not found in mesh!");
386 1686 : }
387 : else
388 : {
389 : // Tally over all mesh blocks if no blocks are provided.
390 2446 : for (const auto & s : _openmc_problem.getMooseMesh().meshSubdomains())
391 1329 : _tally_blocks.insert(s);
392 : }
393 :
394 2803 : _openmc_problem.checkDuplicateEntries(_tally_name, "name");
395 2801 : _openmc_problem.checkDuplicateEntries(_tally_score, "score");
396 :
397 2799 : _local_sum_tally.resize(_tally_score.size(), 0.0);
398 2799 : _local_mean_tally.resize(_tally_score.size(), 0.0);
399 :
400 2799 : _current_tally.resize(_tally_score.size());
401 2799 : _current_raw_tally.resize(_tally_score.size());
402 2799 : _current_raw_tally_rel_error.resize(_tally_score.size());
403 2799 : _current_raw_tally_std_dev.resize(_tally_score.size());
404 2799 : _previous_tally.resize(_tally_score.size());
405 5049 : }
406 :
407 : void
408 3266 : TallyBase::initializeTally()
409 : {
410 : // Clear cached results.
411 3266 : _local_sum_tally.clear();
412 3266 : _local_sum_tally.resize(_tally_score.size(), 0.0);
413 3266 : _local_mean_tally.clear();
414 3266 : _local_mean_tally.resize(_tally_score.size(), 0.0);
415 :
416 3266 : if (_linked_tallies.size() > 0)
417 : {
418 532 : _linked_local_sum_tally.clear();
419 532 : _linked_local_sum_tally.resize(_tally_score.size(), 0.0);
420 : }
421 :
422 3266 : _current_tally.resize(_tally_score.size());
423 3266 : _current_raw_tally.resize(_tally_score.size());
424 3266 : _current_raw_tally_rel_error.resize(_tally_score.size());
425 3266 : _current_raw_tally_std_dev.resize(_tally_score.size());
426 3266 : _previous_tally.resize(_tally_score.size());
427 :
428 3266 : if (_needs_global_tally)
429 : {
430 1658 : _global_sum_tally.clear();
431 1658 : _global_sum_tally.resize(_tally_score.size(), 0.0);
432 : }
433 :
434 : // create the global tally for normalization; we make sure to use the
435 : // same estimator as the local tally
436 3266 : if (addingGlobalTally())
437 : {
438 1528 : _global_tally_index = openmc::model::tallies.size();
439 1528 : _global_tally = openmc::Tally::create();
440 1528 : _global_tally->set_scores(_tally_score);
441 1528 : _global_tally->estimator_ = _estimator;
442 :
443 : // Add a universe filter with the root universe if running in random ray mode.
444 1528 : if (_openmc_problem.runRandomRay())
445 : {
446 : // These vectors are required for OpenMC::span.
447 68 : std::vector<openmc::Filter *> uni_filter = {openmc::Filter::create("universe")};
448 68 : std::vector<int> root_uni = {openmc::model::root_universe};
449 68 : dynamic_cast<openmc::UniverseFilter *>(uni_filter.back())->set_universes(root_uni);
450 68 : _global_tally->set_filters(uni_filter);
451 68 : }
452 : }
453 :
454 3266 : auto [index, spatial_filter] = spatialFilter();
455 3256 : _filter_index = index;
456 :
457 : std::vector<openmc::Filter *> filters;
458 4304 : for (auto & filter : _ext_filters)
459 1048 : filters.push_back(filter->getWrappedFilter());
460 : // We add the spatial filter last to minimize the number of cache
461 : // misses during the OpenMC -> Cardinal transfer.
462 3256 : filters.push_back(spatial_filter);
463 :
464 : // Create the tally, assign the required filters and apply the triggers.
465 3256 : _local_tally_index = openmc::model::tallies.size();
466 3256 : _local_tally = openmc::Tally::create();
467 3256 : _local_tally->set_scores(_tally_score);
468 3256 : _local_tally->estimator_ = _estimator;
469 3256 : _local_tally->set_filters(filters);
470 3256 : applyTriggersToLocalTally(_local_tally);
471 3256 : }
472 :
473 : void
474 624 : TallyBase::resetTally()
475 : {
476 : // Erase the tally.
477 624 : openmc::model::tallies.erase(openmc::model::tallies.begin() + _local_tally_index);
478 :
479 : // Erase global tally.
480 624 : if (addingGlobalTally())
481 150 : openmc::model::tallies.erase(openmc::model::tallies.begin() + _global_tally_index);
482 :
483 : // Erase the filter(s). When running the random ray solver with global normalization we
484 : // add an additional universe filter to the problem before adding the filter at `_filter_index`.
485 : // We need to erase that filter as well so we offset by 1.
486 : const auto erase_idx =
487 624 : addingGlobalTally() && _openmc_problem.runRandomRay() ? _filter_index - 1 : _filter_index;
488 : openmc::model::tally_filters.erase(openmc::model::tally_filters.begin() + erase_idx);
489 624 : }
490 :
491 : Real
492 4680 : TallyBase::storeResults(const std::vector<unsigned int> & var_numbers,
493 : unsigned int local_score,
494 : const std::string & output_type)
495 : {
496 4680 : Real total = 0.0;
497 :
498 4680 : if (output_type == "relaxed")
499 4030 : total += storeResultsInner(var_numbers, local_score, _current_tally);
500 650 : else if (output_type == "rel_error")
501 214 : storeResultsInner(var_numbers, local_score, _current_raw_tally_rel_error, false);
502 436 : else if (output_type == "std_dev")
503 330 : storeResultsInner(var_numbers, local_score, _current_raw_tally_std_dev);
504 106 : else if (output_type == "raw")
505 106 : storeResultsInner(var_numbers, local_score, _current_raw_tally);
506 : else
507 0 : mooseError("Unknown external output " + output_type);
508 :
509 : // Check the normalization.
510 4680 : if (output_type == "relaxed")
511 4030 : checkNormalization(total, local_score);
512 :
513 4680 : return total;
514 : }
515 :
516 : void
517 16 : TallyBase::addScore(const std::string & score)
518 : {
519 16 : _tally_score.push_back(score);
520 :
521 48 : std::vector<std::string> score_names({score});
522 : std::replace(score_names.back().begin(), score_names.back().end(), '-', '_');
523 :
524 : // Modify the variable name and add extra names for the external filter bins.
525 24 : for (const auto & filter : _ext_filters)
526 : {
527 : std::vector<std::string> n;
528 16 : for (unsigned int i = 0; i < score_names.size(); ++i)
529 24 : for (unsigned int j = 0; j < filter->numBins(); ++j)
530 32 : n.push_back(score_names[i] + "_" + filter->binName(j));
531 :
532 8 : score_names = n;
533 8 : }
534 16 : std::copy(score_names.begin(), score_names.end(), std::back_inserter(_tally_name));
535 :
536 16 : _local_sum_tally.resize(_tally_score.size(), 0.0);
537 16 : _local_mean_tally.resize(_tally_score.size(), 0.0);
538 :
539 16 : _current_tally.resize(_tally_score.size());
540 16 : _current_raw_tally.resize(_tally_score.size());
541 16 : _current_raw_tally_rel_error.resize(_tally_score.size());
542 16 : _current_raw_tally_std_dev.resize(_tally_score.size());
543 16 : _previous_tally.resize(_tally_score.size());
544 32 : }
545 :
546 : void
547 2787 : TallyBase::setRelaxation(relaxation::RelaxationEnum relaxation_type, const Real & relaxation_factor)
548 : {
549 2787 : _relaxation_type = relaxation_type;
550 2787 : _relaxation_factor = relaxation_factor;
551 2787 : }
552 :
553 : void
554 3520 : TallyBase::computeSumAndMean()
555 : {
556 7562 : for (unsigned int score = 0; score < _tally_score.size(); ++score)
557 : {
558 4042 : _local_sum_tally[score] = 0.0;
559 :
560 4042 : const unsigned int mapped_bins = _local_tally->n_filter_bins() / _num_ext_filter_bins;
561 9394 : for (unsigned int ext = 0; ext < _num_ext_filter_bins; ++ext)
562 454676 : for (unsigned int m = 0; m < mapped_bins; ++m)
563 449324 : if (!_ext_bins_to_skip[ext])
564 893064 : _local_sum_tally[score] += _local_tally->results_.slice(
565 : openmc::tensor::all,
566 : score,
567 446532 : static_cast<int>(openmc::TallyResult::SUM))[ext * mapped_bins + m];
568 :
569 4042 : _local_mean_tally[score] = _local_sum_tally[score] / _local_tally->n_realizations_;
570 4042 : if (addingGlobalTally())
571 2062 : _global_sum_tally[score] = _openmc_problem.tallySumAcrossBins({_global_tally}, score);
572 :
573 4042 : if (_linked_tallies.size() > 0)
574 864 : _linked_local_sum_tally[score] = _local_sum_tally[score];
575 : }
576 3520 : }
577 :
578 : void
579 2477 : TallyBase::gatherLinkedSum()
580 : {
581 2477 : if (_linked_tallies.size() == 0)
582 : return;
583 :
584 0 : for (const auto & other : _linked_tallies)
585 0 : for (unsigned int score = 0; score < _tally_score.size(); ++score)
586 0 : _linked_local_sum_tally[score] += other->getSum(score);
587 : }
588 :
589 : void
590 3520 : TallyBase::renormalizeLinkedTallies()
591 : {
592 3520 : if (_linked_tallies.size() == 0)
593 : return;
594 :
595 1680 : for (unsigned int score = 0; score < _tally_score.size(); ++score)
596 : {
597 864 : _local_sum_tally[score] = _linked_local_sum_tally[score];
598 864 : _local_mean_tally[score] = _local_sum_tally[score] / _local_tally->n_realizations_;
599 : }
600 : }
601 :
602 : void
603 3518 : TallyBase::relaxAndNormalizeTally()
604 : {
605 : Real alpha;
606 3518 : switch (_relaxation_type)
607 : {
608 : case relaxation::none:
609 : {
610 : alpha = 1.0;
611 : break;
612 : }
613 806 : case relaxation::constant:
614 : {
615 806 : alpha = _relaxation_factor;
616 806 : break;
617 : }
618 24 : case relaxation::robbins_monro:
619 : {
620 24 : alpha = 1.0 / (_openmc_problem.fixedPointIteration() + 1);
621 24 : break;
622 : }
623 48 : case relaxation::dufek_gudowski:
624 : {
625 48 : alpha = static_cast<float>(_openmc_problem.nParticles()) /
626 48 : static_cast<float>(_openmc_problem.nTotalParticles());
627 48 : break;
628 : }
629 0 : default:
630 0 : mooseError("Unhandled RelaxationEnum in TallyBase!");
631 : }
632 :
633 7548 : for (unsigned int score = 0; score < _tally_score.size(); ++score)
634 : {
635 4034 : if (_check_tally_sum && _needs_global_tally)
636 1580 : checkTallySum(score);
637 :
638 4030 : const Real norm = tallyNormalization(score);
639 :
640 4030 : auto & current = _current_tally[score];
641 : auto & previous = _previous_tally[score];
642 : auto & current_raw = _current_raw_tally[score];
643 : auto & current_raw_rel_error = _current_raw_tally_rel_error[score];
644 : auto & current_raw_std_dev = _current_raw_tally_std_dev[score];
645 :
646 4030 : auto mean_tally = _openmc_problem.tallySum(_local_tally, score);
647 : /**
648 : * If the value over the whole domain is zero, then the values in the individual bins must be
649 : * zero. We need to avoid divide-by-zeros.
650 : */
651 4030 : current_raw = mean_tally;
652 4030 : current_raw *= std::abs(norm) < ZERO_TALLY_THRESHOLD ? 0.0 : (1.0 / norm);
653 :
654 4030 : auto sum_sq = OMCTensor(_local_tally->results_.slice(
655 4030 : openmc::tensor::all, score, static_cast<int>(openmc::TallyResult::SUM_SQ)));
656 : current_raw_rel_error =
657 4030 : _openmc_problem.relativeError(mean_tally, sum_sq, _local_tally->n_realizations_);
658 4030 : current_raw_std_dev = current_raw_rel_error * current_raw;
659 :
660 7639 : if (_openmc_problem.fixedPointIteration() == 0 || alpha == 1.0)
661 : {
662 3609 : current = current_raw;
663 3609 : previous = current_raw;
664 : continue;
665 : }
666 :
667 : // Save the current tally (from the previous iteration) into the previous one.
668 : std::copy(current.cbegin(), current.cend(), previous.begin());
669 :
670 : // Relax the tallies by alpha. TODO: skip relaxation when alpha is one.
671 1263 : auto relaxed_tally = (1.0 - alpha) * previous + alpha * current_raw;
672 : std::copy(relaxed_tally.cbegin(), relaxed_tally.cend(), current.begin());
673 : }
674 3514 : }
675 :
676 : void
677 1104 : TallyBase::addLinkedTally(const TallyBase * other)
678 : {
679 1104 : if (this != other)
680 1104 : _linked_tallies.push_back(other);
681 : else
682 0 : mooseError("Internal error: cannot link a tally with itself!");
683 1104 : }
684 :
685 : const openmc::Tally *
686 3636 : TallyBase::getWrappedTally() const
687 : {
688 3636 : if (!_local_tally)
689 0 : mooseError("Internal error: this tally has not been initialized!");
690 :
691 3636 : return _local_tally;
692 : }
693 :
694 : const openmc::Tally *
695 1790 : TallyBase::getWrappedGlobalTally() const
696 : {
697 1790 : if (!_global_tally)
698 0 : mooseError("Internal error: this tally has not been initialized!");
699 :
700 1790 : return _global_tally;
701 : }
702 :
703 : int32_t
704 3254 : TallyBase::getTallyID() const
705 : {
706 3254 : return getWrappedTally()->id();
707 : }
708 :
709 : int32_t
710 1518 : TallyBase::getGlobalTallyID() const
711 : {
712 1518 : return getWrappedGlobalTally()->id();
713 : }
714 :
715 : int
716 1287 : TallyBase::scoreIndex(const std::string & score) const
717 : {
718 1287 : if (!hasScore(score))
719 0 : mooseError("Internal error: tally " + name() + " does not contain the score " + score);
720 :
721 1287 : return std::find(_tally_score.begin(), _tally_score.end(), score) - _tally_score.begin();
722 : }
723 :
724 : std::vector<std::string>
725 152 : TallyBase::getScoreVars(const std::string & score) const
726 : {
727 : std::vector<std::string> score_vars;
728 152 : if (!hasScore(score))
729 : return score_vars;
730 :
731 : unsigned int idx =
732 152 : std::find(_tally_score.begin(), _tally_score.end(), score) - _tally_score.begin();
733 152 : std::copy(_tally_name.begin() + idx * _num_ext_filter_bins,
734 152 : _tally_name.begin() + (idx + 1) * _num_ext_filter_bins,
735 : std::back_inserter(score_vars));
736 :
737 : return score_vars;
738 0 : }
739 :
740 : void
741 489696 : TallyBase::fillElementalAuxVariable(const unsigned int & var_num,
742 : const std::vector<unsigned int> & elem_ids,
743 : const Real & value)
744 : {
745 489696 : auto & solution = _aux.solution();
746 489696 : auto sys_number = _aux.number();
747 :
748 : // loop over all the elements and set the specified variable to the specified value
749 7989046 : for (const auto & e : elem_ids)
750 : {
751 7499350 : auto elem_ptr = _openmc_problem.getMooseMesh().queryElemPtr(e);
752 :
753 7499350 : if (!_openmc_problem.isLocalElem(elem_ptr))
754 3728552 : continue;
755 :
756 3770798 : auto dof_idx = elem_ptr->dof_number(sys_number, var_num, 0);
757 3770798 : solution.set(dof_idx, value);
758 : }
759 489696 : }
760 :
761 : void
762 3256 : TallyBase::applyTriggersToLocalTally(openmc::Tally * tally)
763 : {
764 3256 : if (_tally_trigger)
765 232 : for (int score = 0; score < _tally_score.size(); ++score)
766 256 : tally->triggers_.push_back({_openmc_problem.triggerMetric((*_tally_trigger)[score]),
767 : _tally_trigger_threshold[score],
768 : _trigger_ignore_zeros[score],
769 : score});
770 3256 : }
771 :
772 : Real
773 8060 : TallyBase::tallyNormalization(unsigned int score) const
774 : {
775 8060 : return _normalize_by_global ? _global_sum_tally[score] : _local_sum_tally[score];
776 : }
777 :
778 : void
779 1580 : TallyBase::checkTallySum(const unsigned int & score) const
780 : {
781 : auto rel_diff =
782 1580 : std::abs(_global_sum_tally[score] - _local_sum_tally[score]) / _global_sum_tally[score];
783 1580 : if (rel_diff > _tally_sum_tol)
784 : {
785 4 : std::stringstream msg;
786 4 : bool fixed_source = _openmc_problem.runMode() == openmc::RunMode::FIXED_SOURCE;
787 :
788 8 : msg << _tally_score[score] << " tallies do not match the global " << _tally_score[score]
789 : << " tally:\n"
790 4 : << " Global value: " << Moose::stringify(_global_sum_tally[score])
791 8 : << "\n Tally sum: " << Moose::stringify(_local_sum_tally[score])
792 16 : << "\n Absolute Difference: " << _global_sum_tally[score] - _local_sum_tally[score]
793 4 : << "\n Relative Difference: " << rel_diff;
794 :
795 6 : std::string run_type = fixed_source ? "source_strength" : "power";
796 : msg << "\n\nThis means that the tallies created by Cardinal are missing some hits over the "
797 : "domain. The '"
798 : << run_type
799 : << "' parameter you provided for normalizing the tallies will therefore correspond only to "
800 : "those hits which Cardinal is seeing. If this is your intention, you can set "
801 : "'normalize_by_global_tally' to false, which means that Cardinal will be normalizing "
802 8 : "your tallies so that they correspond to a ";
803 : std::replace(run_type.begin(), run_type.end(), '_', ' ');
804 :
805 : msg << run_type << " of "
806 4 : << Moose::stringify(_local_sum_tally[score] / _global_sum_tally[score] *
807 4 : _openmc_problem.tallyNormalizationValue())
808 12 : << " (NOT the user-specified " << run_type << " of "
809 4 : << _openmc_problem.tallyNormalizationValue()
810 : << "). Or, if this relative difference is acceptable to you, you can loosen "
811 4 : "'tally_sum_tol'.\n\n";
812 :
813 : msg << "If this error is NOT expected, you need to debug your model to determine where some "
814 : "tally hits are not being mapped to Cardinal. If you are using a cell tally, this "
815 : "usually happens because some OpenMC cells are not hit by an element centroid in the "
816 : "[Mesh]. Another possibility is that your OpenMC model is bigger than the [Mesh] and "
817 : "some tally hits are happening outside the [Mesh]. To debug, you can try to:\n\n"
818 : << " 1. Temporarily set 'check_tally_sum = false'\n"
819 : << " 2. Create a MeshTally over the entire domain.\n"
820 : << " 3. Compare the mesh tally to the problematic tally. Do you see any tally hits "
821 4 : "happening in regions where the problematic tally is not picking anything up?";
822 :
823 4 : mooseError(msg.str());
824 0 : }
825 1576 : }
826 :
827 : void
828 4030 : TallyBase::checkNormalization(const Real & sum, unsigned int score) const
829 : {
830 4030 : if (tallyNormalization(score) > ZERO_TALLY_THRESHOLD)
831 4022 : if (_check_tally_sum && std::abs(sum - 1.0) > 1e-6)
832 0 : mooseError("Tally normalization process failed for " + _tally_score[score] +
833 0 : " score! Total fraction of " + Moose::stringify(sum) + " does not match 1.0!");
834 4030 : }
835 : #endif
|
