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/settings.h"
32 :
33 : InputParameters
34 5094 : TallyBase::validParams()
35 : {
36 5094 : auto params = MooseObject::validParams();
37 10188 : params.addParam<MultiMooseEnum>(
38 : "score",
39 10188 : getTallyScoreEnum(),
40 : "Score(s) to use in the OpenMC tallies. If not specified, defaults to 'kappa_fission'");
41 10188 : params.addParam<MooseEnum>(
42 10188 : "estimator", getTallyEstimatorEnum(), "Type of tally estimator to use in OpenMC");
43 10188 : params.addParam<std::vector<std::string>>(
44 : "name",
45 : "Auxiliary variable name(s) to use for OpenMC tallies. "
46 : "If not specified, defaults to the names of the scores");
47 :
48 10188 : params.addParam<std::vector<SubdomainName>>(
49 : "block",
50 : "Subdomains for which to add tallies in OpenMC. If not provided, "
51 : "tallies will be applied over the entire domain corresponding to the [Mesh] block.");
52 10188 : params.addParam<std::vector<SubdomainName>>("blocks",
53 : "This parameter is deprecated, use 'block' instead!");
54 :
55 5094 : MultiMooseEnum tally_trigger("rel_err none");
56 10188 : params.addParam<MultiMooseEnum>(
57 : "trigger",
58 : tally_trigger,
59 : "Trigger criterion to determine when OpenMC simulation is complete "
60 : "based on tallies. If multiple scores are specified in 'score, "
61 : "this same trigger is applied to all scores.");
62 10188 : params.addRangeCheckedParam<std::vector<Real>>(
63 : "trigger_threshold", "trigger_threshold > 0", "Threshold for the tally trigger");
64 10188 : params.addParam<std::vector<bool>>(
65 : "trigger_ignore_zeros",
66 : {false},
67 : "Whether tally bins with zero scores are ignored when computing the tally trigger. If only "
68 : "one "
69 : "value of 'trigger_ignore_zeros' is provided, that value is applied to all tally scores.");
70 :
71 : MultiMooseEnum openmc_outputs(
72 5094 : "unrelaxed_tally_std_dev unrelaxed_tally_rel_error unrelaxed_tally");
73 10188 : params.addParam<MultiMooseEnum>(
74 : "output",
75 : openmc_outputs,
76 : "UNRELAXED field(s) to output from OpenMC for each tally score. "
77 : "unrelaxed_tally_std_dev will write the standard deviation of "
78 : "each tally into auxiliary variables "
79 : "named *_std_dev. unrelaxed_tally_rel_error will write the "
80 : "relative standard deviation (unrelaxed_tally_std_dev / unrelaxed_tally) "
81 : "of each tally into auxiliary variables named *_rel_error. "
82 : "unrelaxed_tally will write the raw unrelaxed tally into auxiliary "
83 : "variables named *_raw (replace * with 'name').");
84 :
85 10188 : params.addParam<std::vector<std::string>>("filters", "External filters to add to this tally.");
86 :
87 10188 : params.addParam<bool>("check_tally_sum",
88 : "Whether to check consistency between the local tallies "
89 : "with a global tally sum. This will require that the "
90 : "integral of the local tally matches a tally with no filters "
91 : "(defined over the entire phase space).");
92 10188 : params.addParam<bool>(
93 : "normalize_by_global_tally",
94 10188 : true,
95 : "Whether to normalize local tallies by a global tally (true) or else by the sum "
96 : "of the local tally (false)");
97 :
98 5094 : params.addPrivateParam<OpenMCCellAverageProblem *>("_openmc_problem");
99 :
100 5094 : params.registerBase("Tally");
101 5094 : params.registerSystemAttributeName("Tally");
102 :
103 5094 : return params;
104 5094 : }
105 :
106 2701 : TallyBase::TallyBase(const InputParameters & parameters)
107 : : MooseObject(parameters),
108 2701 : _openmc_problem(*getParam<OpenMCCellAverageProblem *>("_openmc_problem")),
109 2701 : _mesh(_openmc_problem.mesh()),
110 2701 : _aux(_openmc_problem.getAuxiliarySystem()),
111 2929 : _tally_trigger(isParamValid("trigger") ? &getParam<MultiMooseEnum>("trigger") : nullptr),
112 8103 : _trigger_ignore_zeros(getParam<std::vector<bool>>("trigger_ignore_zeros")),
113 2701 : _normalize_by_global(_openmc_problem.runMode() == openmc::RunMode::FIXED_SOURCE
114 2701 : ? false
115 5296 : : getParam<bool>("normalize_by_global_tally")),
116 2701 : _check_tally_sum(isParamValid("check_tally_sum")
117 6013 : ? getParam<bool>("check_tally_sum")
118 2090 : : (_openmc_problem.runMode() == openmc::RunMode::FIXED_SOURCE
119 2090 : ? true
120 2012 : : _normalize_by_global)),
121 2701 : _needs_global_tally(_check_tally_sum || _normalize_by_global),
122 5402 : _renames_tally_vars(isParamValid("name")),
123 5402 : _has_outputs(isParamValid("output")),
124 8103 : _is_adaptive(_openmc_problem.hasAdaptivity())
125 : {
126 5402 : if (isParamValid("score"))
127 : {
128 3226 : const auto & scores = getParam<MultiMooseEnum>("score");
129 3642 : for (const auto & score : scores)
130 4058 : _tally_score.push_back(_openmc_problem.enumToTallyScore(score));
131 : }
132 : else
133 2176 : _tally_score = {"kappa-fission"};
134 :
135 : const bool heating =
136 2701 : std::find(_tally_score.begin(), _tally_score.end(), "heating") != _tally_score.end();
137 : const bool nu_scatter =
138 2701 : std::find(_tally_score.begin(), _tally_score.end(), "nu-scatter") != _tally_score.end();
139 :
140 5402 : if (isParamValid("estimator"))
141 : {
142 310 : auto estimator = getParam<MooseEnum>("estimator").getEnum<tally::TallyEstimatorEnum>();
143 :
144 : // Photon heating tallies cannot use tracklength estimators.
145 310 : if (estimator == tally::tracklength && openmc::settings::photon_transport && heating)
146 2 : paramError("estimator",
147 : "Tracklength estimators are currently incompatible with photon transport and "
148 : "heating scores! For more information: https://tinyurl.com/3wre3kwt");
149 :
150 308 : if (estimator != tally::analog && nu_scatter)
151 2 : paramError("estimator", "Non-analog estimators are not supported for nu_scatter scores!");
152 :
153 306 : _estimator = _openmc_problem.tallyEstimator(estimator);
154 : }
155 : else
156 : {
157 : /**
158 : * Set a default of tracklength for all tallies other then heating tallies in photon transport
159 : * and nu_scatter tallies. This behavior must be overridden in derived tallies that implement
160 : * mesh filters.
161 : */
162 2391 : _estimator = openmc::TallyEstimator::TRACKLENGTH;
163 :
164 2391 : if (nu_scatter && !(heating && openmc::settings::photon_transport))
165 8 : _estimator = openmc::TallyEstimator::ANALOG;
166 2383 : else if (nu_scatter && heating && openmc::settings::photon_transport)
167 2 : paramError(
168 : "estimator",
169 : "A single tally cannot score both nu_scatter and heating when photon transport is "
170 : "enabled, as both scores require different estimators. Consider adding one tally "
171 : "which scores nu_scatter (with an analog estimator), and a second tally that scores "
172 : "heating (with a collision estimator).");
173 :
174 2389 : if (heating && openmc::settings::photon_transport)
175 8 : _estimator = openmc::TallyEstimator::COLLISION;
176 : }
177 :
178 2695 : if (heating && !openmc::settings::photon_transport)
179 170 : mooseWarning(
180 : "When using the 'heating' score with photon transport disabled, energy deposition\n"
181 : "from photons is neglected unless you specifically ran NJOY to produce MT=301 with\n"
182 : "photon energy deposited locally (not true for any pre-packaged OpenMC data libraries\n"
183 : "on openmc.org).\n\n"
184 : "If you did NOT specifically run NJOY yourself with this customization, we recommend\n"
185 : "using the 'heating_local' score instead, which will capture photon energy deposition.\n"
186 : "Otherwise, you will underpredict the true energy deposition.");
187 :
188 8085 : if (isParamValid("trigger") != isParamValid("trigger_threshold"))
189 2 : paramError("trigger",
190 : "You must either specify none or both of 'trigger' and "
191 : "'trigger_threshold'. You have specified only one.");
192 :
193 2693 : if (_tally_trigger)
194 : {
195 224 : checkRequiredParam(parameters, "trigger_threshold", "using tally triggers");
196 224 : _tally_trigger_threshold = getParam<std::vector<Real>>("trigger_threshold");
197 :
198 112 : if (_tally_trigger->size() != _tally_score.size())
199 2 : paramError("trigger",
200 2 : "'trigger' (size " + std::to_string(_tally_trigger->size()) +
201 2 : ") must have the same length as 'score' (size " +
202 2 : std::to_string(_tally_score.size()) + ")");
203 :
204 110 : if (_tally_trigger_threshold.size() != _tally_score.size())
205 2 : paramError("trigger_threshold",
206 2 : "'trigger_threshold' (size " + std::to_string(_tally_trigger_threshold.size()) +
207 2 : ") must have the same length as 'score' (size " +
208 2 : std::to_string(_tally_score.size()) + ")");
209 :
210 108 : if (_trigger_ignore_zeros.size() > 1)
211 : {
212 2 : if (_tally_score.size() != _trigger_ignore_zeros.size())
213 2 : paramError("trigger_ignore_zeros",
214 2 : "'trigger_ignore_zeros' (size " + std::to_string(_trigger_ignore_zeros.size()) +
215 2 : ") must have the same length as 'score' (size " +
216 2 : std::to_string(_tally_score.size()) + ")");
217 : }
218 106 : else if (_trigger_ignore_zeros.size() == 1)
219 104 : _trigger_ignore_zeros.resize(_tally_score.size(), _trigger_ignore_zeros[0]);
220 :
221 210 : _openmc_problem.checkEmptyVector(_trigger_ignore_zeros, "trigger_ignore_zeros");
222 : }
223 :
224 : // Fetch the filters required by this tally. Error if the filter hasn't been added yet.
225 5370 : if (isParamValid("filters"))
226 : {
227 2292 : for (const auto & filter_name : getParam<std::vector<std::string>>("filters"))
228 : {
229 2042 : if (!_openmc_problem.hasFilter(filter_name))
230 4 : paramError("filters", "Filter with the name " + filter_name + " does not exist!");
231 :
232 1020 : _ext_filters.push_back(_openmc_problem.getFilter(filter_name));
233 : }
234 : }
235 :
236 : // Check the estimator to make sure it doesn't conflict with certain filters.
237 3697 : for (auto & f : _ext_filters)
238 : {
239 1020 : if ((dynamic_cast<AngularLegendreFilter *>(f.get()) ||
240 1020 : dynamic_cast<EnergyOutFilter *>(f.get())) &&
241 186 : _estimator != openmc::TallyEstimator::ANALOG)
242 8 : paramError("estimator",
243 4 : "The filter " + f->name() +
244 : " requires an analog estimator! Please ensure 'estimator' is set to analog.");
245 :
246 1016 : if (dynamic_cast<DelayedGroupFilter *>(f.get()))
247 26 : for (const auto & s : _tally_score)
248 18 : if (s != "delayed-nu-fission" && s != "decay-rate")
249 4 : paramError("score",
250 2 : "The filter " + f->name() +
251 : " can only be used with delayed_nu_fission and decay_rate scores!");
252 : }
253 :
254 5354 : if (isParamValid("name"))
255 1746 : _tally_name = getParam<std::vector<std::string>>("name");
256 : else
257 : {
258 4462 : for (auto score : _tally_score)
259 : {
260 : std::replace(score.begin(), score.end(), '-', '_');
261 2367 : _tally_name.push_back(score);
262 : }
263 : }
264 :
265 2677 : if (_has_outputs)
266 : {
267 : // names of output are appended to ends of 'name'
268 509 : for (const auto & o : getParam<MultiMooseEnum>("output"))
269 : {
270 : std::string name = o;
271 :
272 175 : if (o == "UNRELAXED_TALLY_STD_DEV")
273 132 : _output_name.push_back("std_dev");
274 109 : else if (o == "UNRELAXED_TALLY_REL_ERROR")
275 136 : _output_name.push_back("rel_error");
276 41 : else if (o == "UNRELAXED_TALLY")
277 82 : _output_name.push_back("raw");
278 : else
279 0 : mooseError("Unhandled OutputEnum in OpenMCCellAverageProblem!");
280 : }
281 : }
282 :
283 2677 : if (_tally_name.size() != _tally_score.size())
284 2 : paramError("name", "'name' must be the same length as 'score'!");
285 :
286 : // Modify the variable names so they take into account the bins in the external filters.
287 2675 : auto all_var_names = _tally_name;
288 3689 : for (const auto & filter : _ext_filters)
289 : {
290 : std::vector<std::string> n;
291 2832 : for (unsigned int i = 0; i < all_var_names.size(); ++i)
292 4938 : for (unsigned int j = 0; j < filter->numBins(); ++j)
293 6240 : n.push_back(all_var_names[i] + "_" + filter->binName(j));
294 :
295 1014 : all_var_names = n;
296 :
297 1014 : _num_ext_filter_bins *= filter->numBins();
298 1014 : }
299 2675 : _tally_name = all_var_names;
300 :
301 : // A map of external filter bins to skip when computing sums and means for normalization.
302 2675 : std::vector<bool> skip{false};
303 3689 : for (const auto & filter : _ext_filters)
304 : {
305 : std::vector<bool> s;
306 2610 : for (unsigned int i = 0; i < skip.size(); ++i)
307 4356 : for (unsigned int j = 0; j < filter->numBins(); ++j)
308 3184 : s.push_back(skip[i] || filter->skipBin(j));
309 :
310 1014 : skip = s;
311 : }
312 2675 : _ext_bins_to_skip = skip;
313 :
314 5350 : if (isParamSetByUser("blocks"))
315 0 : mooseError("This parameter is deprecated, use 'block' instead!");
316 :
317 5350 : if (isParamValid("block"))
318 : {
319 4866 : auto block_names = getParam<std::vector<SubdomainName>>("block");
320 1622 : if (block_names.empty())
321 2 : paramError("block", "Subdomain names must be provided if using 'block'!");
322 :
323 1620 : auto block_ids = _openmc_problem.getMooseMesh().getSubdomainIDs(block_names);
324 1620 : std::copy(
325 1620 : block_ids.begin(), block_ids.end(), std::inserter(_tally_blocks, _tally_blocks.end()));
326 :
327 : // Check to make sure all of the blocks are in the mesh.
328 1620 : const auto & subdomains = _openmc_problem.getMooseMesh().meshSubdomains();
329 4133 : for (std::size_t b = 0; b < block_names.size(); ++b)
330 2515 : if (subdomains.find(block_ids[b]) == subdomains.end())
331 4 : paramError("block",
332 2 : "Block '" + block_names[b] + "' specified in 'block' not found in mesh!");
333 1618 : }
334 : else
335 : {
336 : // Tally over all mesh blocks if no blocks are provided.
337 2280 : for (const auto & s : _openmc_problem.getMooseMesh().meshSubdomains())
338 1227 : _tally_blocks.insert(s);
339 : }
340 :
341 2671 : _openmc_problem.checkDuplicateEntries(_tally_name, "name");
342 2669 : _openmc_problem.checkDuplicateEntries(_tally_score, "score");
343 :
344 2667 : _local_sum_tally.resize(_tally_score.size(), 0.0);
345 2667 : _local_mean_tally.resize(_tally_score.size(), 0.0);
346 :
347 2667 : _current_tally.resize(_tally_score.size());
348 2667 : _current_raw_tally.resize(_tally_score.size());
349 2667 : _current_raw_tally_rel_error.resize(_tally_score.size());
350 2667 : _current_raw_tally_std_dev.resize(_tally_score.size());
351 2667 : _previous_tally.resize(_tally_score.size());
352 4843 : }
353 :
354 : void
355 2792 : TallyBase::initializeTally()
356 : {
357 : // Clear cached results.
358 2792 : _local_sum_tally.clear();
359 2792 : _local_sum_tally.resize(_tally_score.size(), 0.0);
360 2792 : _local_mean_tally.clear();
361 2792 : _local_mean_tally.resize(_tally_score.size(), 0.0);
362 :
363 2792 : if (_linked_tallies.size() > 0)
364 : {
365 532 : _linked_local_sum_tally.clear();
366 532 : _linked_local_sum_tally.resize(_tally_score.size(), 0.0);
367 : }
368 :
369 2792 : _current_tally.resize(_tally_score.size());
370 2792 : _current_raw_tally.resize(_tally_score.size());
371 2792 : _current_raw_tally_rel_error.resize(_tally_score.size());
372 2792 : _current_raw_tally_std_dev.resize(_tally_score.size());
373 2792 : _previous_tally.resize(_tally_score.size());
374 :
375 2792 : if (_needs_global_tally)
376 : {
377 1588 : _global_sum_tally.clear();
378 1588 : _global_sum_tally.resize(_tally_score.size(), 0.0);
379 : }
380 :
381 : // create the global tally for normalization; we make sure to use the
382 : // same estimator as the local tally
383 2792 : if (addingGlobalTally())
384 : {
385 1458 : _global_tally_index = openmc::model::tallies.size();
386 1458 : _global_tally = openmc::Tally::create();
387 1458 : _global_tally->set_scores(_tally_score);
388 1458 : _global_tally->estimator_ = _estimator;
389 : }
390 :
391 2792 : auto [index, spatial_filter] = spatialFilter();
392 2782 : _filter_index = index;
393 :
394 : std::vector<openmc::Filter *> filters;
395 3796 : for (auto & filter : _ext_filters)
396 1014 : filters.push_back(filter->getWrappedFilter());
397 : // We add the spatial filter last to minimize the number of cache
398 : // misses during the OpenMC -> Cardinal transfer.
399 2782 : filters.push_back(spatial_filter);
400 :
401 : // Create the tally, assign the required filters and apply the triggers.
402 2782 : _local_tally_index = openmc::model::tallies.size();
403 2782 : _local_tally = openmc::Tally::create();
404 2782 : _local_tally->set_scores(_tally_score);
405 2782 : _local_tally->estimator_ = _estimator;
406 2782 : _local_tally->set_filters(filters);
407 2782 : applyTriggersToLocalTally(_local_tally);
408 2782 : }
409 :
410 : void
411 278 : TallyBase::resetTally()
412 : {
413 : // Erase the tally.
414 278 : openmc::model::tallies.erase(openmc::model::tallies.begin() + _local_tally_index);
415 :
416 : // Erase global tally.
417 278 : if (addingGlobalTally())
418 150 : openmc::model::tallies.erase(openmc::model::tallies.begin() + _global_tally_index);
419 :
420 : // Erase the filter(s).
421 278 : openmc::model::tally_filters.erase(openmc::model::tally_filters.begin() + _filter_index);
422 278 : }
423 :
424 : Real
425 4066 : TallyBase::storeResults(const std::vector<unsigned int> & var_numbers,
426 : unsigned int local_score,
427 : const std::string & output_type)
428 : {
429 4066 : Real total = 0.0;
430 :
431 4066 : if (output_type == "relaxed")
432 3770 : total += storeResultsInner(var_numbers, local_score, _current_tally);
433 296 : else if (output_type == "rel_error")
434 78 : storeResultsInner(var_numbers, local_score, _current_raw_tally_rel_error, false);
435 218 : else if (output_type == "std_dev")
436 112 : storeResultsInner(var_numbers, local_score, _current_raw_tally_std_dev);
437 106 : else if (output_type == "raw")
438 106 : storeResultsInner(var_numbers, local_score, _current_raw_tally);
439 : else
440 0 : mooseError("Unknown external output " + output_type);
441 :
442 : // Check the normalization.
443 4066 : if (output_type == "relaxed")
444 3770 : checkNormalization(total, local_score);
445 :
446 4066 : return total;
447 : }
448 :
449 : void
450 16 : TallyBase::addScore(const std::string & score)
451 : {
452 16 : _tally_score.push_back(score);
453 :
454 48 : std::vector<std::string> score_names({score});
455 : std::replace(score_names.back().begin(), score_names.back().end(), '-', '_');
456 :
457 : // Modify the variable name and add extra names for the external filter bins.
458 24 : for (const auto & filter : _ext_filters)
459 : {
460 : std::vector<std::string> n;
461 16 : for (unsigned int i = 0; i < score_names.size(); ++i)
462 24 : for (unsigned int j = 0; j < filter->numBins(); ++j)
463 32 : n.push_back(score_names[i] + "_" + filter->binName(j));
464 :
465 8 : score_names = n;
466 8 : }
467 16 : std::copy(score_names.begin(), score_names.end(), std::back_inserter(_tally_name));
468 :
469 16 : _local_sum_tally.resize(_tally_score.size(), 0.0);
470 16 : _local_mean_tally.resize(_tally_score.size(), 0.0);
471 :
472 16 : _current_tally.resize(_tally_score.size());
473 16 : _current_raw_tally.resize(_tally_score.size());
474 16 : _current_raw_tally_rel_error.resize(_tally_score.size());
475 16 : _current_raw_tally_std_dev.resize(_tally_score.size());
476 16 : _previous_tally.resize(_tally_score.size());
477 32 : }
478 :
479 : void
480 2657 : TallyBase::setRelaxation(relaxation::RelaxationEnum relaxation_type, const Real & relaxation_factor)
481 : {
482 2657 : _relaxation_type = relaxation_type;
483 2657 : _relaxation_factor = relaxation_factor;
484 2657 : }
485 :
486 : void
487 3328 : TallyBase::computeSumAndMean()
488 : {
489 7100 : for (unsigned int score = 0; score < _tally_score.size(); ++score)
490 : {
491 3772 : _local_sum_tally[score] = 0.0;
492 :
493 3772 : const unsigned int mapped_bins = _local_tally->n_filter_bins() / _num_ext_filter_bins;
494 8854 : for (unsigned int ext = 0; ext < _num_ext_filter_bins; ++ext)
495 452402 : for (unsigned int m = 0; m < mapped_bins; ++m)
496 447320 : if (!_ext_bins_to_skip[ext])
497 889056 : _local_sum_tally[score] += _local_tally->results_.slice(
498 : openmc::tensor::all,
499 : score,
500 444528 : static_cast<int>(openmc::TallyResult::SUM))[ext * mapped_bins + m];
501 :
502 3772 : _local_mean_tally[score] = _local_sum_tally[score] / _local_tally->n_realizations_;
503 3772 : if (addingGlobalTally())
504 1914 : _global_sum_tally[score] = _openmc_problem.tallySumAcrossBins({_global_tally}, score);
505 :
506 3772 : if (_linked_tallies.size() > 0)
507 864 : _linked_local_sum_tally[score] = _local_sum_tally[score];
508 : }
509 3328 : }
510 :
511 : void
512 2287 : TallyBase::gatherLinkedSum()
513 : {
514 2287 : if (_linked_tallies.size() == 0)
515 : return;
516 :
517 0 : for (const auto & other : _linked_tallies)
518 0 : for (unsigned int score = 0; score < _tally_score.size(); ++score)
519 0 : _linked_local_sum_tally[score] += other->getSum(score);
520 : }
521 :
522 : void
523 3328 : TallyBase::renormalizeLinkedTallies()
524 : {
525 3328 : if (_linked_tallies.size() == 0)
526 : return;
527 :
528 1680 : for (unsigned int score = 0; score < _tally_score.size(); ++score)
529 : {
530 864 : _local_sum_tally[score] = _linked_local_sum_tally[score];
531 864 : _local_mean_tally[score] = _local_sum_tally[score] / _local_tally->n_realizations_;
532 : }
533 : }
534 :
535 : void
536 3328 : TallyBase::relaxAndNormalizeTally()
537 : {
538 : Real alpha;
539 3328 : switch (_relaxation_type)
540 : {
541 : case relaxation::none:
542 : {
543 : alpha = 1.0;
544 : break;
545 : }
546 734 : case relaxation::constant:
547 : {
548 734 : alpha = _relaxation_factor;
549 734 : break;
550 : }
551 24 : case relaxation::robbins_monro:
552 : {
553 24 : alpha = 1.0 / (_openmc_problem.fixedPointIteration() + 1);
554 24 : break;
555 : }
556 48 : case relaxation::dufek_gudowski:
557 : {
558 48 : alpha = static_cast<float>(_openmc_problem.nParticles()) /
559 48 : static_cast<float>(_openmc_problem.nTotalParticles());
560 48 : break;
561 : }
562 0 : default:
563 0 : mooseError("Unhandled RelaxationEnum in TallyBase!");
564 : }
565 :
566 7098 : for (unsigned int score = 0; score < _tally_score.size(); ++score)
567 : {
568 3772 : if (_check_tally_sum && _needs_global_tally)
569 1442 : checkTallySum(score);
570 :
571 3770 : const Real norm = tallyNormalization(score);
572 :
573 3770 : auto & current = _current_tally[score];
574 : auto & previous = _previous_tally[score];
575 : auto & current_raw = _current_raw_tally[score];
576 : auto & current_raw_rel_error = _current_raw_tally_rel_error[score];
577 : auto & current_raw_std_dev = _current_raw_tally_std_dev[score];
578 :
579 3770 : auto mean_tally = _openmc_problem.tallySum(_local_tally, score);
580 : /**
581 : * If the value over the whole domain is zero, then the values in the individual bins must be
582 : * zero. We need to avoid divide-by-zeros.
583 : */
584 3770 : current_raw = mean_tally;
585 3770 : current_raw *= std::abs(norm) < ZERO_TALLY_THRESHOLD ? 0.0 : (1.0 / norm);
586 :
587 3770 : auto sum_sq = OMCTensor(_local_tally->results_.slice(
588 3770 : openmc::tensor::all, score, static_cast<int>(openmc::TallyResult::SUM_SQ)));
589 : current_raw_rel_error =
590 3770 : _openmc_problem.relativeError(mean_tally, sum_sq, _local_tally->n_realizations_);
591 3770 : current_raw_std_dev = current_raw_rel_error * current_raw;
592 :
593 7167 : if (_openmc_problem.fixedPointIteration() == 0 || alpha == 1.0)
594 : {
595 3397 : current = current_raw;
596 3397 : previous = current_raw;
597 : continue;
598 : }
599 :
600 : // Save the current tally (from the previous iteration) into the previous one.
601 : std::copy(current.cbegin(), current.cend(), previous.begin());
602 :
603 : // Relax the tallies by alpha. TODO: skip relaxation when alpha is one.
604 1119 : auto relaxed_tally = (1.0 - alpha) * previous + alpha * current_raw;
605 : std::copy(relaxed_tally.cbegin(), relaxed_tally.cend(), current.begin());
606 : }
607 3326 : }
608 :
609 : void
610 1104 : TallyBase::addLinkedTally(const TallyBase * other)
611 : {
612 1104 : if (this != other)
613 1104 : _linked_tallies.push_back(other);
614 : else
615 0 : mooseError("Internal error: cannot link a tally with itself!");
616 1104 : }
617 :
618 : const openmc::Tally *
619 3160 : TallyBase::getWrappedTally() const
620 : {
621 3160 : if (!_local_tally)
622 0 : mooseError("Internal error: this tally has not been initialized!");
623 :
624 3160 : return _local_tally;
625 : }
626 :
627 : const openmc::Tally *
628 1720 : TallyBase::getWrappedGlobalTally() const
629 : {
630 1720 : if (!_global_tally)
631 0 : mooseError("Internal error: this tally has not been initialized!");
632 :
633 1720 : return _global_tally;
634 : }
635 :
636 : int32_t
637 2780 : TallyBase::getTallyID() const
638 : {
639 2780 : return getWrappedTally()->id();
640 : }
641 :
642 : int32_t
643 1448 : TallyBase::getGlobalTallyID() const
644 : {
645 1448 : return getWrappedGlobalTally()->id();
646 : }
647 :
648 : int
649 1256 : TallyBase::scoreIndex(const std::string & score) const
650 : {
651 1256 : if (!hasScore(score))
652 0 : mooseError("Internal error: tally " + name() + " does not contain the score " + score);
653 :
654 1256 : return std::find(_tally_score.begin(), _tally_score.end(), score) - _tally_score.begin();
655 : }
656 :
657 : std::vector<std::string>
658 152 : TallyBase::getScoreVars(const std::string & score) const
659 : {
660 : std::vector<std::string> score_vars;
661 152 : if (!hasScore(score))
662 : return score_vars;
663 :
664 : unsigned int idx =
665 152 : std::find(_tally_score.begin(), _tally_score.end(), score) - _tally_score.begin();
666 152 : std::copy(_tally_name.begin() + idx * _num_ext_filter_bins,
667 152 : _tally_name.begin() + (idx + 1) * _num_ext_filter_bins,
668 : std::back_inserter(score_vars));
669 :
670 : return score_vars;
671 0 : }
672 :
673 : void
674 485486 : TallyBase::fillElementalAuxVariable(const unsigned int & var_num,
675 : const std::vector<unsigned int> & elem_ids,
676 : const Real & value)
677 : {
678 485486 : auto & solution = _aux.solution();
679 485486 : auto sys_number = _aux.number();
680 :
681 : // loop over all the elements and set the specified variable to the specified value
682 7961342 : for (const auto & e : elem_ids)
683 : {
684 7475856 : auto elem_ptr = _openmc_problem.getMooseMesh().queryElemPtr(e);
685 :
686 7475856 : if (!_openmc_problem.isLocalElem(elem_ptr))
687 3722360 : continue;
688 :
689 3753496 : auto dof_idx = elem_ptr->dof_number(sys_number, var_num, 0);
690 3753496 : solution.set(dof_idx, value);
691 : }
692 485486 : }
693 :
694 : void
695 2782 : TallyBase::applyTriggersToLocalTally(openmc::Tally * tally)
696 : {
697 2782 : if (_tally_trigger)
698 232 : for (int score = 0; score < _tally_score.size(); ++score)
699 256 : tally->triggers_.push_back({_openmc_problem.triggerMetric((*_tally_trigger)[score]),
700 : _tally_trigger_threshold[score],
701 : _trigger_ignore_zeros[score],
702 : score});
703 2782 : }
704 :
705 : Real
706 7540 : TallyBase::tallyNormalization(unsigned int score) const
707 : {
708 7540 : return _normalize_by_global ? _global_sum_tally[score] : _local_sum_tally[score];
709 : }
710 :
711 : void
712 1442 : TallyBase::checkTallySum(const unsigned int & score) const
713 : {
714 1442 : if (std::abs(_global_sum_tally[score] - _local_sum_tally[score]) / _global_sum_tally[score] >
715 : 1e-6)
716 : {
717 2 : std::stringstream msg;
718 4 : msg << _tally_score[score] << " tallies do not match the global " << _tally_score[score]
719 : << " tally:\n"
720 2 : << " Global value: " << Moose::stringify(_global_sum_tally[score])
721 4 : << "\n Tally sum: " << Moose::stringify(_local_sum_tally[score])
722 8 : << "\n Difference: " << _global_sum_tally[score] - _local_sum_tally[score]
723 : << "\n\nThis means that the tallies created by Cardinal are missing some hits over the "
724 : "domain.\n"
725 2 : << "You can turn off this check by setting 'check_tally_sum' to false.";
726 :
727 2 : mooseError(msg.str());
728 0 : }
729 1440 : }
730 :
731 : void
732 3770 : TallyBase::checkNormalization(const Real & sum, unsigned int score) const
733 : {
734 3770 : if (tallyNormalization(score) > ZERO_TALLY_THRESHOLD)
735 3762 : if (_check_tally_sum && std::abs(sum - 1.0) > 1e-6)
736 0 : mooseError("Tally normalization process failed for " + _tally_score[score] +
737 0 : " score! Total fraction of " + Moose::stringify(sum) + " does not match 1.0!");
738 3770 : }
739 : #endif
|
