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 "MeshTally.h"
21 : #include "DisplacedProblem.h"
22 :
23 : #include "libmesh/replicated_mesh.h"
24 :
25 : registerMooseObject("CardinalApp", MeshTally);
26 :
27 : InputParameters
28 1161 : MeshTally::validParams()
29 : {
30 1161 : auto params = TallyBase::validParams();
31 1161 : params.addClassDescription("A class which implements unstructured mesh tallies.");
32 2322 : params.addParam<std::string>("mesh_template",
33 : "Mesh tally template for OpenMC when using mesh tallies; "
34 : "at present, this mesh must exactly match the mesh used in the "
35 : "[Mesh] block because a one-to-one copy "
36 : "is used to get OpenMC's tally results on the [Mesh].");
37 2322 : params.addParam<Point>("mesh_translation",
38 : "Coordinate to which this mesh should be "
39 : "translated. Units must match those used to define the [Mesh].");
40 :
41 : // The index of this tally into an array of mesh translations. Defaults to zero.
42 1161 : params.addPrivateParam<unsigned int>("instance", 0);
43 :
44 1161 : return params;
45 0 : }
46 :
47 767 : MeshTally::MeshTally(const InputParameters & parameters)
48 : : TallyBase(parameters),
49 2082 : _mesh_translation(isParamValid("mesh_translation") ? getParam<Point>("mesh_translation")
50 : : Point(0.0, 0.0, 0.0)),
51 1530 : _instance(getParam<unsigned int>("instance")),
52 2974 : _use_dof_map(_is_adaptive || isParamValid("block"))
53 : {
54 : bool nu_scatter =
55 765 : std::find(_tally_score.begin(), _tally_score.end(), "nu-scatter") != _tally_score.end();
56 :
57 : // Error check the estimators.
58 1530 : if (isParamValid("estimator"))
59 : {
60 14 : if (_estimator == openmc::TallyEstimator::TRACKLENGTH)
61 2 : paramError("estimator",
62 : "Tracklength estimators are currently incompatible with mesh tallies!");
63 : }
64 : else
65 1502 : _estimator = nu_scatter ? openmc::TallyEstimator::ANALOG : openmc::TallyEstimator::COLLISION;
66 :
67 : // Error check the mesh template.
68 765 : if (_openmc_problem.getMooseMesh().getMesh().allow_renumbering() &&
69 2 : !_openmc_problem.getMooseMesh().getMesh().is_replicated())
70 2 : mooseError(
71 : "Mesh tallies currently require 'allow_renumbering = false' to be set in the [Mesh]!");
72 :
73 1522 : if (isParamValid("mesh_template"))
74 : {
75 652 : if (_is_adaptive)
76 2 : paramError("mesh_template",
77 : "Adaptivity is not supported when loading a mesh from 'mesh_template'!");
78 :
79 1300 : if (isParamValid("block"))
80 2 : paramError("block",
81 : "Block restriction is currently not supported for mesh tallies which load a "
82 : "mesh from a file!");
83 :
84 648 : if (_openmc_problem.useDisplaced())
85 1 : paramError("mesh_template",
86 : "Tallying on a file-based mesh is not supported for moving-mesh cases as there is "
87 : "not a mechanism to update the mesh geometry. You can use a mesh tally for moving "
88 : "mesh cases by instead tallying directly on the [Mesh]. Simply do not provide the "
89 : "'mesh_template' parameter.");
90 :
91 1294 : _mesh_template_filename = &getParam<std::string>("mesh_template");
92 : }
93 : else
94 : {
95 : // for distributed meshes, each rank only owns a portion of the mesh information, but
96 : // OpenMC wants the entire mesh to be available on every rank. We might be able to add
97 : // this feature in the future, but will need to investigate
98 109 : if (!_openmc_problem.getMooseMesh().getMesh().is_replicated())
99 0 : mooseError("Directly tallying on the [Mesh] block by OpenMC is not yet supported "
100 : "for distributed meshes!");
101 :
102 218 : if (isParamValid("mesh_translation"))
103 0 : paramError("mesh_translation",
104 : "The mesh filter cannot be translated if directly tallying on the mesh "
105 : "provided in the [Mesh] block!");
106 : }
107 :
108 : /**
109 : * If the instance isn't zero this variable is a translated mesh tally. It will accumulate it's
110 : * scores in a different set of variables (the auxvars which are added by the first tally in a
111 : * sequence of mesh tallies), and so it doesn't need to create any auxvars.
112 : */
113 756 : if (_instance != 0)
114 368 : _tally_name = std::vector<std::string>();
115 756 : }
116 :
117 : std::pair<unsigned int, openmc::Filter *>
118 769 : MeshTally::spatialFilter()
119 : {
120 : // Create the OpenMC mesh which will be tallied on.
121 769 : if (!_mesh_template_filename)
122 : {
123 : auto msh =
124 137 : dynamic_cast<const libMesh::ReplicatedMesh *>(_openmc_problem.getMooseMesh().getMeshPtr());
125 137 : if (!msh)
126 0 : mooseError("Internal error: The mesh is not a replicated mesh.");
127 :
128 : // Adaptivity and block restriction both require a map from the element subsets we want to
129 : // tally on to the full mesh.
130 137 : if (_use_dof_map)
131 : {
132 86 : _bin_to_element_mapping.clear();
133 :
134 : auto begin = _tally_blocks.size() > 0
135 86 : ? msh->active_subdomain_set_elements_begin(_tally_blocks)
136 172 : : msh->active_elements_begin();
137 86 : auto end = _tally_blocks.size() > 0 ? msh->active_subdomain_set_elements_end(_tally_blocks)
138 86 : : msh->active_elements_end();
139 25772 : for (const auto & old_elem : libMesh::as_range(begin, end))
140 25686 : _bin_to_element_mapping.push_back(old_elem->id());
141 :
142 : _bin_to_element_mapping.shrink_to_fit();
143 : }
144 :
145 137 : openmc::model::meshes.emplace_back(std::make_unique<openmc::AdaptiveLibMesh>(
146 137 : _openmc_problem.getMooseMesh().getMesh(), _openmc_problem.scaling(), _tally_blocks));
147 : }
148 : else
149 632 : openmc::model::meshes.emplace_back(
150 1264 : std::make_unique<openmc::LibMesh>(*_mesh_template_filename, _openmc_problem.scaling()));
151 :
152 769 : _mesh_index = openmc::model::meshes.size() - 1;
153 769 : _mesh_template =
154 769 : dynamic_cast<openmc::UnstructuredMesh *>(openmc::model::meshes[_mesh_index].get());
155 :
156 : // by setting the ID to -1, OpenMC will automatically detect the next available ID
157 769 : _mesh_template->set_id(-1);
158 769 : _mesh_template->output_ = false;
159 :
160 769 : _mesh_filter = dynamic_cast<openmc::MeshFilter *>(openmc::Filter::create("mesh"));
161 769 : _mesh_filter->set_mesh(_mesh_index);
162 769 : _mesh_filter->set_translation({_mesh_translation(0), _mesh_translation(1), _mesh_translation(2)});
163 :
164 : // Validate the mesh filters to make sure we can run a copy transfer to the [Mesh].
165 769 : checkMeshTemplateAndTranslations();
166 :
167 763 : return std::make_pair(openmc::model::tally_filters.size() - 1, _mesh_filter);
168 : }
169 :
170 : void
171 61 : MeshTally::resetTally()
172 : {
173 61 : TallyBase::resetTally();
174 :
175 : // Erase the OpenMC mesh.
176 61 : openmc::model::meshes.erase(openmc::model::meshes.begin() + _mesh_index);
177 61 : }
178 :
179 : void
180 1041 : MeshTally::gatherLinkedSum()
181 : {
182 1041 : if (_linked_tallies.size() == 0)
183 : return;
184 :
185 2448 : for (const auto & other : _linked_tallies)
186 : {
187 3360 : for (unsigned int score = 0; score < _tally_score.size(); ++score)
188 : {
189 1728 : _linked_local_sum_tally[score] += other->getSum(score);
190 1728 : if (other->addingGlobalTally())
191 272 : _global_sum_tally[score] =
192 272 : _openmc_problem.tallySumAcrossBins({other->getWrappedGlobalTally()}, score);
193 : }
194 : }
195 : }
196 :
197 : Real
198 1253 : MeshTally::storeResultsInner(const std::vector<unsigned int> & var_numbers,
199 : unsigned int local_score,
200 : const std::vector<OMCTensor> & tally_vals,
201 : bool norm_by_src_rate)
202 : {
203 : Real total = 0.0;
204 :
205 1253 : unsigned int mesh_offset = _instance * _mesh_filter->n_bins();
206 2752 : for (unsigned int ext_bin = 0; ext_bin < _num_ext_filter_bins; ++ext_bin)
207 : {
208 469979 : for (decltype(_mesh_filter->n_bins()) e = 0; e < _mesh_filter->n_bins(); ++e)
209 : {
210 468480 : Real unnormalized_tally = tally_vals[local_score](ext_bin * _mesh_filter->n_bins() + e);
211 :
212 : // divide each tally by the volume that it corresponds to in MOOSE
213 : // because we will apply it as a volumetric tally (per unit volume).
214 : // Because we require that the mesh template has units of cm based on the
215 : // mesh constructors in OpenMC, we need to adjust the division
216 468480 : Real volumetric_tally = unnormalized_tally;
217 468480 : volumetric_tally *= norm_by_src_rate
218 923332 : ? _openmc_problem.tallyMultiplier(_tally_score[local_score],
219 454852 : _local_mean_tally[local_score]) /
220 454852 : _mesh_template->volume(e) * _openmc_problem.scaling() *
221 : _openmc_problem.scaling() * _openmc_problem.scaling()
222 : : 1.0;
223 468480 : total += _ext_bins_to_skip[ext_bin] ? 0.0 : unnormalized_tally;
224 :
225 468480 : auto var = var_numbers[_num_ext_filter_bins * local_score + ext_bin];
226 468480 : auto elem_id = _use_dof_map ? _bin_to_element_mapping[e] : mesh_offset + e;
227 468480 : fillElementalAuxVariable(var, {elem_id}, volumetric_tally);
228 : }
229 : }
230 :
231 1253 : return total;
232 : }
233 :
234 : void
235 769 : MeshTally::checkMeshTemplateAndTranslations()
236 : {
237 : // we can do some rudimentary checking on the mesh template by comparing the centroid
238 : // coordinates compared to centroids in the [Mesh] (because right now, we just doing a simple
239 : // copy transfer that necessitates the meshes to have the same elements in the same order). In
240 : // other words, you might have two meshes that represent the same geometry, the element ordering
241 : // could be different.
242 769 : unsigned int mesh_offset = _instance * _mesh_filter->n_bins();
243 242361 : for (int e = 0; e < _mesh_filter->n_bins(); ++e)
244 : {
245 241598 : auto elem_id = _use_dof_map ? _bin_to_element_mapping[e] : mesh_offset + e;
246 241598 : auto elem_ptr = _openmc_problem.getMooseMesh().queryElemPtr(elem_id);
247 :
248 : // if element is not on this part of the distributed mesh, skip it
249 241598 : if (!elem_ptr)
250 51402 : continue;
251 :
252 190196 : const auto pt = _mesh_template->centroid(e);
253 190196 : Point centroid_template = {pt[0], pt[1], pt[2]};
254 :
255 : // The translation applied in OpenMC isn't actually registered in the mesh itself;
256 : // it is always added on to the point, so we need to do the same here
257 : centroid_template += _mesh_translation;
258 :
259 : // because the mesh template and [Mesh] may be in different units, we need
260 : // to adjust the [Mesh] by the scaling factor before doing a comparison.
261 190196 : Point centroid_mesh = elem_ptr->vertex_average() * _openmc_problem.scaling();
262 :
263 : // if the centroids are the same except for a factor of 'scaling', then we can
264 : // guess that the mesh_template is probably not in units of centimeters
265 190196 : if (_openmc_problem.hasScaling())
266 : {
267 : // if scaling was applied correctly, then each calculation of 'scaling' here should equal 1.
268 : // Otherwise, if they're all the same, then 'scaling_x' is probably the factor by which the
269 : // mesh_template needs to be multiplied, so we can print a helpful error message
270 : bool incorrect_scaling = true;
271 122024 : for (unsigned int j = 0; j < OpenMCCellAverageProblem::DIMENSION; ++j)
272 : {
273 91518 : Real scaling = centroid_mesh(j) / centroid_template(j);
274 91518 : incorrect_scaling = incorrect_scaling && !MooseUtils::absoluteFuzzyEqual(scaling, 1.0);
275 : }
276 :
277 30506 : if (incorrect_scaling)
278 2 : paramError("mesh_template",
279 : "The centroids of the 'mesh_template' differ from the "
280 2 : "centroids of the [Mesh] by a factor of " +
281 2 : Moose::stringify(centroid_mesh(0) / centroid_template(0)) +
282 : ".\nDid you forget that the 'mesh_template' must be in "
283 : "the same units as the [Mesh]?");
284 : }
285 :
286 : // check if centroids are the same
287 : bool different_centroids = false;
288 760776 : for (unsigned int j = 0; j < OpenMCCellAverageProblem::DIMENSION; ++j)
289 570582 : different_centroids = different_centroids ||
290 : !MooseUtils::absoluteFuzzyEqual(centroid_mesh(j), centroid_template(j));
291 :
292 190194 : if (different_centroids)
293 4 : paramError(
294 : "mesh_template",
295 4 : "Centroid for element " + Moose::stringify(elem_id) +
296 8 : " in the [Mesh] (cm): " + _openmc_problem.printPoint(centroid_mesh) +
297 4 : "\ndoes not match centroid for element " + Moose::stringify(e) +
298 8 : " in the 'mesh_template' with instance " + Moose::stringify(_instance) +
299 8 : " (cm): " + _openmc_problem.printPoint(centroid_template) +
300 : "!\n\nThe copy transfer requires that the [Mesh] and 'mesh_template' be identical.");
301 : }
302 763 : }
303 : #endif
|
