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 :
30 : #include "openmc/settings.h"
31 :
32 : InputParameters
33 4785 : TallyBase::validParams()
34 : {
35 4785 : auto params = MooseObject::validParams();
36 9570 : params.addParam<MultiMooseEnum>(
37 : "score",
38 9570 : getTallyScoreEnum(),
39 : "Score(s) to use in the OpenMC tallies. If not specified, defaults to 'kappa_fission'");
40 9570 : params.addParam<MooseEnum>(
41 9570 : "estimator", getTallyEstimatorEnum(), "Type of tally estimator to use in OpenMC");
42 9570 : params.addParam<std::vector<std::string>>(
43 : "name",
44 : "Auxiliary variable name(s) to use for OpenMC tallies. "
45 : "If not specified, defaults to the names of the scores");
46 :
47 4785 : MultiMooseEnum tally_trigger("rel_err none");
48 9570 : params.addParam<MultiMooseEnum>(
49 : "trigger",
50 : tally_trigger,
51 : "Trigger criterion to determine when OpenMC simulation is complete "
52 : "based on tallies. If multiple scores are specified in 'score, "
53 : "this same trigger is applied to all scores.");
54 9570 : params.addRangeCheckedParam<std::vector<Real>>(
55 : "trigger_threshold", "trigger_threshold > 0", "Threshold for the tally trigger");
56 9570 : params.addParam<std::vector<bool>>(
57 : "trigger_ignore_zeros",
58 : {false},
59 : "Whether tally bins with zero scores are ignored when computing the tally trigger. If only "
60 : "one "
61 : "value of 'trigger_ignore_zeros' is provided, that value is applied to all tally scores.");
62 :
63 : MultiMooseEnum openmc_outputs(
64 4785 : "unrelaxed_tally_std_dev unrelaxed_tally_rel_error unrelaxed_tally");
65 9570 : params.addParam<MultiMooseEnum>(
66 : "output",
67 : openmc_outputs,
68 : "UNRELAXED field(s) to output from OpenMC for each tally score. "
69 : "unrelaxed_tally_std_dev will write the standard deviation of "
70 : "each tally into auxiliary variables "
71 : "named *_std_dev. unrelaxed_tally_rel_error will write the "
72 : "relative standard deviation (unrelaxed_tally_std_dev / unrelaxed_tally) "
73 : "of each tally into auxiliary variables named *_rel_error. "
74 : "unrelaxed_tally will write the raw unrelaxed tally into auxiliary "
75 : "variables named *_raw (replace * with 'name').");
76 :
77 9570 : params.addParam<std::vector<std::string>>("filters", "External filters to add to this tally.");
78 :
79 4785 : params.addPrivateParam<OpenMCCellAverageProblem *>("_openmc_problem");
80 :
81 4785 : params.registerBase("Tally");
82 4785 : params.registerSystemAttributeName("Tally");
83 :
84 4785 : return params;
85 4785 : }
86 :
87 2540 : TallyBase::TallyBase(const InputParameters & parameters)
88 : : MooseObject(parameters),
89 2540 : _openmc_problem(*getParam<OpenMCCellAverageProblem *>("_openmc_problem")),
90 2540 : _mesh(_openmc_problem.mesh()),
91 2540 : _aux(_openmc_problem.getAuxiliarySystem()),
92 5328 : _tally_trigger(isParamValid("trigger") ? &getParam<MultiMooseEnum>("trigger") : nullptr),
93 5080 : _trigger_ignore_zeros(getParam<std::vector<bool>>("trigger_ignore_zeros")),
94 5080 : _renames_tally_vars(isParamValid("name")),
95 5080 : _has_outputs(isParamValid("output")),
96 12700 : _is_adaptive(_openmc_problem.hasAdaptivity())
97 : {
98 5080 : if (isParamValid("score"))
99 : {
100 2904 : const auto & scores = getParam<MultiMooseEnum>("score");
101 3274 : for (const auto & score : scores)
102 3644 : _tally_score.push_back(_openmc_problem.enumToTallyScore(score));
103 : }
104 : else
105 2176 : _tally_score = {"kappa-fission"};
106 :
107 : bool heating =
108 2540 : std::find(_tally_score.begin(), _tally_score.end(), "heating") != _tally_score.end();
109 : bool nu_scatter =
110 2540 : std::find(_tally_score.begin(), _tally_score.end(), "nu-scatter") != _tally_score.end();
111 :
112 5080 : if (isParamValid("estimator"))
113 : {
114 350 : auto estimator = getParam<MooseEnum>("estimator").getEnum<tally::TallyEstimatorEnum>();
115 :
116 : // Photon heating tallies cannot use tracklength estimators.
117 350 : if (estimator == tally::tracklength && openmc::settings::photon_transport && heating)
118 2 : paramError("estimator",
119 : "Tracklength estimators are currently incompatible with photon transport and "
120 : "heating scores! For more information: https://tinyurl.com/3wre3kwt");
121 :
122 348 : if (estimator != tally::analog && nu_scatter)
123 2 : paramError("estimator", "Non-analog estimators are not supported for nu_scatter scores!");
124 :
125 346 : _estimator = _openmc_problem.tallyEstimator(estimator);
126 : }
127 : else
128 : {
129 : /**
130 : * Set a default of tracklength for all tallies other then heating tallies in photon transport
131 : * and nu_scatter tallies. This behavior must be overridden in derived tallies that implement
132 : * mesh filters.
133 : */
134 2190 : _estimator = openmc::TallyEstimator::TRACKLENGTH;
135 :
136 2190 : if (nu_scatter && !(heating && openmc::settings::photon_transport))
137 8 : _estimator = openmc::TallyEstimator::ANALOG;
138 2182 : else if (nu_scatter && heating && openmc::settings::photon_transport)
139 2 : paramError(
140 : "estimator",
141 : "A single tally cannot score both nu_scatter and heating when photon transport is "
142 : "enabled, as both scores require different estimators. Consider adding one tally "
143 : "which scores nu_scatter (with an analog estimator), and a second tally that scores "
144 : "heating (with a collision estimator).");
145 :
146 2188 : if (heating && openmc::settings::photon_transport)
147 8 : _estimator = openmc::TallyEstimator::COLLISION;
148 : }
149 :
150 2534 : if (heating && !openmc::settings::photon_transport)
151 140 : mooseWarning(
152 : "When using the 'heating' score with photon transport disabled, energy deposition\n"
153 : "from photons is neglected unless you specifically ran NJOY to produce MT=301 with\n"
154 : "photon energy deposited locally (not true for any pre-packaged OpenMC data libraries\n"
155 : "on openmc.org).\n\n"
156 : "If you did NOT specifically run NJOY yourself with this customization, we recommend\n"
157 : "using the 'heating_local' score instead, which will capture photon energy deposition.\n"
158 : "Otherwise, you will underpredict the true energy deposition.");
159 :
160 7602 : if (isParamValid("trigger") != isParamValid("trigger_threshold"))
161 2 : paramError("trigger",
162 : "You must either specify none or both of 'trigger' and "
163 : "'trigger_threshold'. You have specified only one.");
164 :
165 2532 : if (_tally_trigger)
166 : {
167 244 : checkRequiredParam(parameters, "trigger_threshold", "using tally triggers");
168 244 : _tally_trigger_threshold = getParam<std::vector<Real>>("trigger_threshold");
169 :
170 122 : if (_tally_trigger->size() != _tally_score.size())
171 2 : paramError("trigger",
172 2 : "'trigger' (size " + std::to_string(_tally_trigger->size()) +
173 2 : ") must have the same length as 'score' (size " +
174 2 : std::to_string(_tally_score.size()) + ")");
175 :
176 120 : if (_tally_trigger_threshold.size() != _tally_score.size())
177 2 : paramError("trigger_threshold",
178 2 : "'trigger_threshold' (size " + std::to_string(_tally_trigger_threshold.size()) +
179 2 : ") must have the same length as 'score' (size " +
180 2 : std::to_string(_tally_score.size()) + ")");
181 :
182 118 : if (_trigger_ignore_zeros.size() > 1)
183 : {
184 2 : if (_tally_score.size() != _trigger_ignore_zeros.size())
185 2 : paramError("trigger_ignore_zeros",
186 2 : "'trigger_ignore_zeros' (size " + std::to_string(_trigger_ignore_zeros.size()) +
187 2 : ") must have the same length as 'score' (size " +
188 2 : std::to_string(_tally_score.size()) + ")");
189 : }
190 116 : else if (_trigger_ignore_zeros.size() == 1)
191 114 : _trigger_ignore_zeros.resize(_tally_score.size(), _trigger_ignore_zeros[0]);
192 :
193 230 : _openmc_problem.checkEmptyVector(_trigger_ignore_zeros, "trigger_ignore_zeros");
194 : }
195 :
196 : // Fetch the filters required by this tally. Error if the filter hasn't been added yet.
197 5048 : if (isParamValid("filters"))
198 : {
199 2452 : for (const auto & filter_name : getParam<std::vector<std::string>>("filters"))
200 : {
201 2242 : if (!_openmc_problem.hasFilter(filter_name))
202 4 : paramError("filters", "Filter with the name " + filter_name + " does not exist!");
203 :
204 1120 : _ext_filters.push_back(_openmc_problem.getFilter(filter_name));
205 : }
206 : }
207 :
208 : // Check the estimator to make sure it doesn't conflict with certain filters.
209 3638 : for (auto & f : _ext_filters)
210 1120 : if ((dynamic_cast<AngularLegendreFilter *>(f.get()) ||
211 1120 : dynamic_cast<EnergyOutFilter *>(f.get())) &&
212 216 : _estimator != openmc::TallyEstimator::ANALOG)
213 4 : paramError("estimator",
214 4 : "The filter " + f->name() +
215 : " requires an analog estimator! Please ensure 'estimator' is set to analog.");
216 :
217 5036 : if (isParamValid("name"))
218 1716 : _tally_name = getParam<std::vector<std::string>>("name");
219 : else
220 : {
221 4142 : for (auto score : _tally_score)
222 : {
223 : std::replace(score.begin(), score.end(), '-', '_');
224 2196 : _tally_name.push_back(score);
225 : }
226 : }
227 :
228 2518 : if (_has_outputs)
229 : {
230 : // names of output are appended to ends of 'name'
231 422 : for (const auto & o : getParam<MultiMooseEnum>("output"))
232 : {
233 : std::string name = o;
234 :
235 150 : if (o == "UNRELAXED_TALLY_STD_DEV")
236 172 : _output_name.push_back("std_dev");
237 64 : else if (o == "UNRELAXED_TALLY_REL_ERROR")
238 48 : _output_name.push_back("rel_error");
239 40 : else if (o == "UNRELAXED_TALLY")
240 80 : _output_name.push_back("raw");
241 : else
242 0 : mooseError("Unhandled OutputEnum in OpenMCCellAverageProblem!");
243 : }
244 : }
245 :
246 2518 : if (_tally_name.size() != _tally_score.size())
247 2 : paramError("name", "'name' must be the same length as 'score'!");
248 :
249 : // Modify the variable names so they take into account the bins in the external filters.
250 2516 : auto all_var_names = _tally_name;
251 3632 : for (const auto & filter : _ext_filters)
252 : {
253 : std::vector<std::string> n;
254 3188 : for (unsigned int i = 0; i < all_var_names.size(); ++i)
255 5456 : for (unsigned int j = 0; j < filter->numBins(); ++j)
256 6768 : n.push_back(all_var_names[i] + "_" + filter->binName(j));
257 :
258 1116 : all_var_names = n;
259 :
260 1116 : _num_ext_filter_bins *= filter->numBins();
261 1116 : }
262 2516 : _tally_name = all_var_names;
263 :
264 : // A map of external filter bins to skip when computing sums and means for normalization.
265 2516 : std::vector<bool> skip{false};
266 3632 : for (const auto & filter : _ext_filters)
267 : {
268 : std::vector<bool> s;
269 2944 : for (unsigned int i = 0; i < skip.size(); ++i)
270 4860 : for (unsigned int j = 0; j < filter->numBins(); ++j)
271 3456 : s.push_back(skip[i] || filter->skipBin(j));
272 :
273 1116 : skip = s;
274 : }
275 2516 : _ext_bins_to_skip = skip;
276 :
277 2516 : _openmc_problem.checkDuplicateEntries(_tally_name, "name");
278 2514 : _openmc_problem.checkDuplicateEntries(_tally_score, "score");
279 :
280 2512 : _local_sum_tally.resize(_tally_score.size(), 0.0);
281 2512 : _local_mean_tally.resize(_tally_score.size(), 0.0);
282 :
283 2512 : _current_tally.resize(_tally_score.size());
284 2512 : _current_raw_tally.resize(_tally_score.size());
285 2512 : _current_raw_tally_rel_error.resize(_tally_score.size());
286 2512 : _current_raw_tally_std_dev.resize(_tally_score.size());
287 2512 : _previous_tally.resize(_tally_score.size());
288 4688 : }
289 :
290 : void
291 2486 : TallyBase::initializeTally()
292 : {
293 : // Clear cached results.
294 2486 : _local_sum_tally.clear();
295 2486 : _local_sum_tally.resize(_tally_score.size(), 0.0);
296 2486 : _local_mean_tally.clear();
297 2486 : _local_mean_tally.resize(_tally_score.size(), 0.0);
298 :
299 2486 : _current_tally.resize(_tally_score.size());
300 2486 : _current_raw_tally.resize(_tally_score.size());
301 2486 : _current_raw_tally_rel_error.resize(_tally_score.size());
302 2486 : _current_raw_tally_std_dev.resize(_tally_score.size());
303 2486 : _previous_tally.resize(_tally_score.size());
304 :
305 2486 : auto [index, spatial_filter] = spatialFilter();
306 2476 : _filter_index = index;
307 :
308 : std::vector<openmc::Filter *> filters;
309 3592 : for (auto & filter : _ext_filters)
310 1116 : filters.push_back(filter->getWrappedFilter());
311 : // We add the spatial filter last to minimize the number of cache
312 : // misses during the OpenMC -> Cardinal transfer.
313 2476 : filters.push_back(spatial_filter);
314 :
315 : // Create the tally, assign the required filters and apply the triggers.
316 2476 : _local_tally_index = openmc::model::tallies.size();
317 2476 : _local_tally = openmc::Tally::create();
318 2476 : _local_tally->set_scores(_tally_score);
319 2476 : _local_tally->estimator_ = _estimator;
320 2476 : _local_tally->set_filters(filters);
321 2476 : applyTriggersToLocalTally(_local_tally);
322 2476 : }
323 :
324 : void
325 84 : TallyBase::resetTally()
326 : {
327 : // Erase the tally.
328 84 : openmc::model::tallies.erase(openmc::model::tallies.begin() + _local_tally_index);
329 :
330 : // Erase the filter(s).
331 84 : openmc::model::tally_filters.erase(openmc::model::tally_filters.begin() + _filter_index);
332 84 : }
333 :
334 : Real
335 3820 : TallyBase::storeResults(const std::vector<unsigned int> & var_numbers,
336 : unsigned int local_score,
337 : unsigned int global_score,
338 : const std::string & output_type)
339 : {
340 : Real total = 0.0;
341 :
342 3820 : if (output_type == "relaxed")
343 3548 : total += storeResultsInner(var_numbers, local_score, global_score, _current_tally);
344 272 : else if (output_type == "rel_error")
345 32 : storeResultsInner(var_numbers, local_score, global_score, _current_raw_tally_rel_error, false);
346 240 : else if (output_type == "std_dev")
347 136 : storeResultsInner(var_numbers, local_score, global_score, _current_raw_tally_std_dev);
348 104 : else if (output_type == "raw")
349 104 : storeResultsInner(var_numbers, local_score, global_score, _current_raw_tally);
350 : else
351 0 : mooseError("Unknown external output " + output_type);
352 :
353 3820 : return total;
354 : }
355 :
356 : void
357 16 : TallyBase::addScore(const std::string & score)
358 : {
359 16 : _tally_score.push_back(score);
360 :
361 48 : std::vector<std::string> score_names({score});
362 : std::replace(score_names.back().begin(), score_names.back().end(), '-', '_');
363 :
364 : // Modify the variable name and add extra names for the external filter bins.
365 24 : for (const auto & filter : _ext_filters)
366 : {
367 : std::vector<std::string> n;
368 16 : for (unsigned int i = 0; i < score_names.size(); ++i)
369 24 : for (unsigned int j = 0; j < filter->numBins(); ++j)
370 32 : n.push_back(score_names[i] + "_" + filter->binName(j));
371 :
372 8 : score_names = n;
373 8 : }
374 16 : std::copy(score_names.begin(), score_names.end(), std::back_inserter(_tally_name));
375 :
376 16 : _local_sum_tally.resize(_tally_score.size(), 0.0);
377 16 : _local_mean_tally.resize(_tally_score.size(), 0.0);
378 :
379 16 : _current_tally.resize(_tally_score.size());
380 16 : _current_raw_tally.resize(_tally_score.size());
381 16 : _current_raw_tally_rel_error.resize(_tally_score.size());
382 16 : _current_raw_tally_std_dev.resize(_tally_score.size());
383 16 : _previous_tally.resize(_tally_score.size());
384 32 : }
385 :
386 : void
387 3180 : TallyBase::computeSumAndMean()
388 : {
389 6730 : for (unsigned int score = 0; score < _tally_score.size(); ++score)
390 : {
391 3550 : _local_sum_tally[score] = 0.0;
392 :
393 3550 : const unsigned int mapped_bins = _local_tally->n_filter_bins() / _num_ext_filter_bins;
394 8420 : for (unsigned int ext = 0; ext < _num_ext_filter_bins; ++ext)
395 870480 : for (unsigned int m = 0; m < mapped_bins; ++m)
396 865610 : if (!_ext_bins_to_skip[ext])
397 862818 : _local_sum_tally[score] +=
398 1725636 : xt::view(_local_tally->results_,
399 862818 : xt::all(),
400 : score,
401 862818 : static_cast<int>(openmc::TallyResult::SUM))[ext * mapped_bins + m];
402 :
403 3550 : _local_mean_tally[score] = _local_sum_tally[score] / _local_tally->n_realizations_;
404 : }
405 3180 : }
406 :
407 : void
408 3548 : TallyBase::relaxAndNormalizeTally(unsigned int local_score, const Real & alpha, const Real & norm)
409 : {
410 3548 : auto & current = _current_tally[local_score];
411 : auto & previous = _previous_tally[local_score];
412 : auto & current_raw = _current_raw_tally[local_score];
413 : auto & current_raw_rel_error = _current_raw_tally_rel_error[local_score];
414 : auto & current_raw_std_dev = _current_raw_tally_std_dev[local_score];
415 :
416 3548 : auto mean_tally = _openmc_problem.tallySum(_local_tally, local_score);
417 : /**
418 : * If the value over the whole domain is zero, then the values in the individual bins must be
419 : * zero. We need to avoid divide-by-zeros.
420 : */
421 3548 : current_raw = std::abs(norm) < ZERO_TALLY_THRESHOLD
422 10652 : ? static_cast<xt::xtensor<double, 1>>(mean_tally * 0.0)
423 7096 : : static_cast<xt::xtensor<double, 1>>(mean_tally / norm);
424 :
425 3548 : auto sum_sq = xt::view(_local_tally->results_,
426 3548 : xt::all(),
427 : local_score,
428 3548 : static_cast<int>(openmc::TallyResult::SUM_SQ));
429 : current_raw_rel_error =
430 3548 : _openmc_problem.relativeError(mean_tally, sum_sq, _local_tally->n_realizations_);
431 0 : current_raw_std_dev = current_raw_rel_error * current_raw;
432 :
433 3548 : if (_openmc_problem.fixedPointIteration() == 0 || alpha == 1.0)
434 : {
435 3204 : current = current_raw;
436 3204 : previous = current_raw;
437 : return;
438 : }
439 :
440 : // Save the current tally (from the previous iteration) into the previous one.
441 344 : std::copy(current.cbegin(), current.cend(), previous.begin());
442 :
443 : // Relax the tallies by alpha. TODO: skip relaxation when alpha is one.
444 688 : auto relaxed_tally = (1.0 - alpha) * previous + alpha * current_raw;
445 344 : std::copy(relaxed_tally.cbegin(), relaxed_tally.cend(), current.begin());
446 3892 : }
447 :
448 : const openmc::Tally *
449 2976 : TallyBase::getWrappedTally() const
450 : {
451 2976 : if (!_local_tally)
452 0 : mooseError("This tally has not been initialized!");
453 :
454 2976 : return _local_tally;
455 : }
456 :
457 : int32_t
458 2474 : TallyBase::getTallyID() const
459 : {
460 2474 : return getWrappedTally()->id();
461 : }
462 :
463 : int
464 1004 : TallyBase::scoreIndex(const std::string & score) const
465 : {
466 1004 : if (!hasScore(score))
467 0 : mooseError("Internal error: tally " + name() + " does not contain the score " + score);
468 :
469 1004 : return std::find(_tally_score.begin(), _tally_score.end(), score) - _tally_score.begin();
470 : }
471 :
472 : std::vector<std::string>
473 104 : TallyBase::getScoreVars(const std::string & score) const
474 : {
475 : std::vector<std::string> score_vars;
476 104 : if (!hasScore(score))
477 : return score_vars;
478 :
479 : unsigned int idx =
480 104 : std::find(_tally_score.begin(), _tally_score.end(), score) - _tally_score.begin();
481 104 : std::copy(_tally_name.begin() + idx * _num_ext_filter_bins,
482 104 : _tally_name.begin() + (idx + 1) * _num_ext_filter_bins,
483 : std::back_inserter(score_vars));
484 :
485 : return score_vars;
486 0 : }
487 :
488 : void
489 1416594 : TallyBase::fillElementalAuxVariable(const unsigned int & var_num,
490 : const std::vector<unsigned int> & elem_ids,
491 : const Real & value)
492 : {
493 1416594 : auto & solution = _aux.solution();
494 1416594 : auto sys_number = _aux.number();
495 :
496 : // loop over all the elements and set the specified variable to the specified value
497 95966258 : for (const auto & e : elem_ids)
498 : {
499 94549664 : auto elem_ptr = _openmc_problem.getMooseMesh().queryElemPtr(e);
500 :
501 94549664 : if (!_openmc_problem.isLocalElem(elem_ptr))
502 55931984 : continue;
503 :
504 38617680 : auto dof_idx = elem_ptr->dof_number(sys_number, var_num, 0);
505 38617680 : solution.set(dof_idx, value);
506 : }
507 1416594 : }
508 :
509 : void
510 2476 : TallyBase::applyTriggersToLocalTally(openmc::Tally * tally)
511 : {
512 2476 : if (_tally_trigger)
513 252 : for (int score = 0; score < _tally_score.size(); ++score)
514 276 : tally->triggers_.push_back({_openmc_problem.triggerMetric((*_tally_trigger)[score]),
515 : _tally_trigger_threshold[score],
516 : _trigger_ignore_zeros[score],
517 : score});
518 2476 : }
519 : #endif
|
