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 1167 : MeshTally::validParams()
29 : {
30 1167 : auto params = TallyBase::validParams();
31 1167 : params.addClassDescription("A class which implements unstructured mesh tallies.");
32 2334 : 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 2334 : 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 1167 : params.addPrivateParam<unsigned int>("instance", 0);
43 :
44 1167 : return params;
45 0 : }
46 :
47 771 : MeshTally::MeshTally(const InputParameters & parameters)
48 : : TallyBase(parameters),
49 2090 : _mesh_translation(isParamValid("mesh_translation") ? getParam<Point>("mesh_translation")
50 : : Point(0.0, 0.0, 0.0)),
51 1538 : _instance(getParam<unsigned int>("instance")),
52 2990 : _use_dof_map(_is_adaptive || isParamValid("block"))
53 : {
54 : bool nu_scatter =
55 769 : std::find(_tally_score.begin(), _tally_score.end(), "nu-scatter") != _tally_score.end();
56 :
57 : // Error check the estimators.
58 1538 : 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 1510 : _estimator = nu_scatter ? openmc::TallyEstimator::ANALOG : openmc::TallyEstimator::COLLISION;
66 :
67 : // Error check the mesh template.
68 769 : 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 1530 : 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 113 : 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 226 : 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 760 : if (_instance != 0)
114 368 : _tally_name = std::vector<std::string>();
115 :
116 : // The random ray solver requires tracklength estimators, which unstructured meshes don't support.
117 760 : if (_openmc_problem.runRandomRay())
118 2 : mooseError("Unstructured mesh tallies are not supported when using the random ray solver!");
119 758 : }
120 :
121 : std::pair<unsigned int, openmc::Filter *>
122 771 : MeshTally::spatialFilter()
123 : {
124 : // Create the OpenMC mesh which will be tallied on.
125 771 : if (!_mesh_template_filename)
126 : {
127 : auto msh =
128 139 : dynamic_cast<const libMesh::ReplicatedMesh *>(_openmc_problem.getMooseMesh().getMeshPtr());
129 139 : if (!msh)
130 0 : mooseError("Internal error: The mesh is not a replicated mesh.");
131 :
132 : // Adaptivity and block restriction both require a map from the element subsets we want to
133 : // tally on to the full mesh.
134 139 : if (_use_dof_map)
135 : {
136 86 : _bin_to_element_mapping.clear();
137 :
138 : auto begin = _tally_blocks.size() > 0
139 86 : ? msh->active_subdomain_set_elements_begin(_tally_blocks)
140 172 : : msh->active_elements_begin();
141 86 : auto end = _tally_blocks.size() > 0 ? msh->active_subdomain_set_elements_end(_tally_blocks)
142 86 : : msh->active_elements_end();
143 25772 : for (const auto & old_elem : libMesh::as_range(begin, end))
144 25686 : _bin_to_element_mapping.push_back(old_elem->id());
145 :
146 : _bin_to_element_mapping.shrink_to_fit();
147 : }
148 :
149 139 : openmc::model::meshes.emplace_back(std::make_unique<openmc::AdaptiveLibMesh>(
150 139 : _openmc_problem.getMooseMesh().getMesh(), _openmc_problem.scaling(), _tally_blocks));
151 : }
152 : else
153 632 : openmc::model::meshes.emplace_back(
154 1264 : std::make_unique<openmc::LibMesh>(*_mesh_template_filename, _openmc_problem.scaling()));
155 :
156 771 : _mesh_index = openmc::model::meshes.size() - 1;
157 771 : _mesh_template =
158 771 : dynamic_cast<openmc::UnstructuredMesh *>(openmc::model::meshes[_mesh_index].get());
159 :
160 : // by setting the ID to -1, OpenMC will automatically detect the next available ID
161 771 : _mesh_template->set_id(-1);
162 771 : _mesh_template->output_ = false;
163 :
164 771 : _mesh_filter = dynamic_cast<openmc::MeshFilter *>(openmc::Filter::create("mesh"));
165 771 : _mesh_filter->set_mesh(_mesh_index);
166 771 : _mesh_filter->set_translation({_mesh_translation(0), _mesh_translation(1), _mesh_translation(2)});
167 :
168 : // Validate the mesh filters to make sure we can run a copy transfer to the [Mesh].
169 771 : checkMeshTemplateAndTranslations();
170 :
171 765 : return std::make_pair(openmc::model::tally_filters.size() - 1, _mesh_filter);
172 : }
173 :
174 : void
175 61 : MeshTally::resetTally()
176 : {
177 61 : TallyBase::resetTally();
178 :
179 : // Erase the OpenMC mesh.
180 61 : openmc::model::meshes.erase(openmc::model::meshes.begin() + _mesh_index);
181 61 : }
182 :
183 : void
184 1043 : MeshTally::gatherLinkedSum()
185 : {
186 1043 : if (_linked_tallies.size() == 0)
187 : return;
188 :
189 2448 : for (const auto & other : _linked_tallies)
190 : {
191 3360 : for (unsigned int score = 0; score < _tally_score.size(); ++score)
192 : {
193 1728 : _linked_local_sum_tally[score] += other->getSum(score);
194 1728 : if (other->addingGlobalTally())
195 272 : _global_sum_tally[score] =
196 272 : _openmc_problem.tallySumAcrossBins({other->getWrappedGlobalTally()}, score);
197 : }
198 : }
199 : }
200 :
201 : Real
202 1255 : MeshTally::storeResultsInner(const std::vector<unsigned int> & var_numbers,
203 : unsigned int local_score,
204 : const std::vector<OMCTensor> & tally_vals,
205 : bool norm_by_src_rate)
206 : {
207 : Real total = 0.0;
208 :
209 1255 : unsigned int mesh_offset = _instance * _mesh_filter->n_bins();
210 2756 : for (unsigned int ext_bin = 0; ext_bin < _num_ext_filter_bins; ++ext_bin)
211 : {
212 469983 : for (decltype(_mesh_filter->n_bins()) e = 0; e < _mesh_filter->n_bins(); ++e)
213 : {
214 468482 : Real unnormalized_tally = tally_vals[local_score](ext_bin * _mesh_filter->n_bins() + e);
215 :
216 : // divide each tally by the volume that it corresponds to in MOOSE
217 : // because we will apply it as a volumetric tally (per unit volume).
218 : // Because we require that the mesh template has units of cm based on the
219 : // mesh constructors in OpenMC, we need to adjust the division
220 468482 : Real volumetric_tally = unnormalized_tally;
221 468482 : volumetric_tally *= norm_by_src_rate
222 923336 : ? _openmc_problem.tallyMultiplier(_tally_score[local_score],
223 454854 : _local_mean_tally[local_score]) /
224 454854 : _mesh_template->volume(e) * _openmc_problem.scaling() *
225 : _openmc_problem.scaling() * _openmc_problem.scaling()
226 : : 1.0;
227 468482 : total += _ext_bins_to_skip[ext_bin] ? 0.0 : unnormalized_tally;
228 :
229 468482 : auto var = var_numbers[_num_ext_filter_bins * local_score + ext_bin];
230 468482 : auto elem_id = _use_dof_map ? _bin_to_element_mapping[e] : mesh_offset + e;
231 468482 : fillElementalAuxVariable(var, {elem_id}, volumetric_tally);
232 : }
233 : }
234 :
235 1255 : return total;
236 : }
237 :
238 : void
239 771 : MeshTally::checkMeshTemplateAndTranslations()
240 : {
241 : // we can do some rudimentary checking on the mesh template by comparing the centroid
242 : // coordinates compared to centroids in the [Mesh] (because right now, we just doing a simple
243 : // copy transfer that necessitates the meshes to have the same elements in the same order). In
244 : // other words, you might have two meshes that represent the same geometry, the element ordering
245 : // could be different.
246 771 : unsigned int mesh_offset = _instance * _mesh_filter->n_bins();
247 242365 : for (int e = 0; e < _mesh_filter->n_bins(); ++e)
248 : {
249 241600 : auto elem_id = _use_dof_map ? _bin_to_element_mapping[e] : mesh_offset + e;
250 241600 : auto elem_ptr = _openmc_problem.getMooseMesh().queryElemPtr(elem_id);
251 :
252 : // if element is not on this part of the distributed mesh, skip it
253 241600 : if (!elem_ptr)
254 51402 : continue;
255 :
256 190198 : const auto pt = _mesh_template->centroid(e);
257 190198 : Point centroid_template = {pt[0], pt[1], pt[2]};
258 :
259 : // The translation applied in OpenMC isn't actually registered in the mesh itself;
260 : // it is always added on to the point, so we need to do the same here
261 : centroid_template += _mesh_translation;
262 :
263 : // because the mesh template and [Mesh] may be in different units, we need
264 : // to adjust the [Mesh] by the scaling factor before doing a comparison.
265 190198 : Point centroid_mesh = elem_ptr->vertex_average() * _openmc_problem.scaling();
266 :
267 : // if the centroids are the same except for a factor of 'scaling', then we can
268 : // guess that the mesh_template is probably not in units of centimeters
269 190198 : if (_openmc_problem.hasScaling())
270 : {
271 : // if scaling was applied correctly, then each calculation of 'scaling' here should equal 1.
272 : // Otherwise, if they're all the same, then 'scaling_x' is probably the factor by which the
273 : // mesh_template needs to be multiplied, so we can print a helpful error message
274 : bool incorrect_scaling = true;
275 122024 : for (unsigned int j = 0; j < OpenMCCellAverageProblem::DIMENSION; ++j)
276 : {
277 91518 : Real scaling = centroid_mesh(j) / centroid_template(j);
278 91518 : incorrect_scaling = incorrect_scaling && !MooseUtils::absoluteFuzzyEqual(scaling, 1.0);
279 : }
280 :
281 30506 : if (incorrect_scaling)
282 2 : paramError("mesh_template",
283 : "The centroids of the 'mesh_template' differ from the "
284 2 : "centroids of the [Mesh] by a factor of " +
285 2 : Moose::stringify(centroid_mesh(0) / centroid_template(0)) +
286 : ".\nDid you forget that the 'mesh_template' must be in "
287 : "the same units as the [Mesh]?");
288 : }
289 :
290 : // check if centroids are the same
291 : bool different_centroids = false;
292 760784 : for (unsigned int j = 0; j < OpenMCCellAverageProblem::DIMENSION; ++j)
293 570588 : different_centroids = different_centroids ||
294 : !MooseUtils::absoluteFuzzyEqual(centroid_mesh(j), centroid_template(j));
295 :
296 190196 : if (different_centroids)
297 4 : paramError(
298 : "mesh_template",
299 4 : "Centroid for element " + Moose::stringify(elem_id) +
300 8 : " in the [Mesh] (cm): " + _openmc_problem.printPoint(centroid_mesh) +
301 4 : "\ndoes not match centroid for element " + Moose::stringify(e) +
302 8 : " in the 'mesh_template' with instance " + Moose::stringify(_instance) +
303 8 : " (cm): " + _openmc_problem.printPoint(centroid_template) +
304 : "!\n\nThe copy transfer requires that the [Mesh] and 'mesh_template' be identical.");
305 : }
306 765 : }
307 : #endif
|
