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 1089 : MeshTally::validParams()
29 : {
30 1089 : auto params = TallyBase::validParams();
31 1089 : params.addClassDescription("A class which implements unstructured mesh tallies.");
32 2178 : 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 2178 : 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 2178 : params.addParam<std::vector<SubdomainName>>(
41 : "block",
42 : "Subdomains for which to add tallies in OpenMC. If not provided, this mesh "
43 : "tally will be applied over the entire mesh.");
44 2178 : params.addParam<std::vector<SubdomainName>>("blocks",
45 : "This parameter is deprecated, use 'block' instead!");
46 :
47 : // The index of this tally into an array of mesh translations. Defaults to zero.
48 1089 : params.addPrivateParam<unsigned int>("instance", 0);
49 :
50 1089 : return params;
51 0 : }
52 :
53 726 : MeshTally::MeshTally(const InputParameters & parameters)
54 : : TallyBase(parameters),
55 1986 : _mesh_translation(isParamValid("mesh_translation") ? getParam<Point>("mesh_translation")
56 : : Point(0.0, 0.0, 0.0)),
57 1452 : _instance(getParam<unsigned int>("instance")),
58 2836 : _use_dof_map(_is_adaptive || isParamValid("block"))
59 : {
60 1452 : if (isParamSetByUser("blocks"))
61 0 : mooseError("This parameter is deprecated, use 'block' instead!");
62 :
63 : bool nu_scatter =
64 726 : std::find(_tally_score.begin(), _tally_score.end(), "nu-scatter") != _tally_score.end();
65 :
66 : // Error check the estimators.
67 1452 : if (isParamValid("estimator"))
68 : {
69 14 : if (_estimator == openmc::TallyEstimator::TRACKLENGTH)
70 2 : paramError("estimator",
71 : "Tracklength estimators are currently incompatible with mesh tallies!");
72 : }
73 : else
74 1424 : _estimator = nu_scatter ? openmc::TallyEstimator::ANALOG : openmc::TallyEstimator::COLLISION;
75 :
76 : // Error check the mesh template.
77 726 : if (_openmc_problem.getMooseMesh().getMesh().allow_renumbering() &&
78 2 : !_openmc_problem.getMooseMesh().getMesh().is_replicated())
79 2 : mooseError(
80 : "Mesh tallies currently require 'allow_renumbering = false' to be set in the [Mesh]!");
81 :
82 1444 : if (isParamValid("mesh_template"))
83 : {
84 641 : if (_is_adaptive)
85 2 : paramError("mesh_template",
86 : "Adaptivity is not supported when loading a mesh from 'mesh_template'!");
87 :
88 1278 : if (isParamValid("block"))
89 2 : paramError("block",
90 : "Block restriction is currently not supported for mesh tallies which load a "
91 : "mesh from a file!");
92 :
93 637 : if (_openmc_problem.useDisplaced())
94 1 : paramError("mesh_template",
95 : "Tallying on a file-based mesh is not supported for moving-mesh cases as there is "
96 : "not a mechanism to update the mesh geometry. You can use a mesh tally for moving "
97 : "mesh cases by instead tallying directly on the [Mesh]. Simply do not provide the "
98 : "'mesh_template' parameter.");
99 :
100 1272 : _mesh_template_filename = &getParam<std::string>("mesh_template");
101 : }
102 : else
103 : {
104 81 : if (std::abs(_openmc_problem.scaling() - 1.0) > 1e-6)
105 2 : mooseError("Directly tallying on the [Mesh] is only supported for 'scaling' of unity. "
106 : "Instead, please make a file containing your tally mesh and set it with "
107 : "'mesh_template'. You can generate a mesh file corresponding to the [Mesh] "
108 2 : "by running:\n\ncardinal-opt -i " +
109 2 : _app.getFileName() + " --mesh-only");
110 :
111 : // for distributed meshes, each rank only owns a portion of the mesh information, but
112 : // OpenMC wants the entire mesh to be available on every rank. We might be able to add
113 : // this feature in the future, but will need to investigate
114 79 : if (!_openmc_problem.getMooseMesh().getMesh().is_replicated())
115 0 : mooseError("Directly tallying on the [Mesh] block by OpenMC is not yet supported "
116 : "for distributed meshes!");
117 :
118 158 : if (isParamValid("mesh_translation"))
119 0 : paramError("mesh_translation",
120 : "The mesh filter cannot be translated if directly tallying on the mesh "
121 : "provided in the [Mesh] block!");
122 :
123 : // Fetch subdomain IDs for block restrictions.
124 158 : if (isParamValid("block"))
125 : {
126 48 : auto block_names = getParam<std::vector<SubdomainName>>("block");
127 16 : if (block_names.empty())
128 2 : paramError("block", "Subdomain names must be provided if using 'block'!");
129 :
130 14 : auto block_ids = _openmc_problem.getMooseMesh().getSubdomainIDs(block_names);
131 14 : std::copy(
132 : block_ids.begin(), block_ids.end(), std::inserter(_tally_blocks, _tally_blocks.end()));
133 :
134 : // Check to make sure all of the blocks are in the mesh.
135 14 : const auto & subdomains = _openmc_problem.getMooseMesh().meshSubdomains();
136 34 : for (std::size_t b = 0; b < block_names.size(); ++b)
137 20 : if (subdomains.find(block_ids[b]) == subdomains.end())
138 0 : paramError("block",
139 0 : "Block '" + block_names[b] + "' specified in 'block' not found in mesh!");
140 14 : }
141 : }
142 :
143 : /**
144 : * If the instance isn't zero this variable is a translated mesh tally. It will accumulate it's
145 : * scores in a different set of variables (the auxvars which are added by the first tally in a
146 : * sequence of mesh tallies), and so it doesn't need to create any auxvars.
147 : */
148 713 : if (_instance != 0)
149 356 : _tally_name = std::vector<std::string>();
150 713 : }
151 :
152 : std::pair<unsigned int, openmc::Filter *>
153 722 : MeshTally::spatialFilter()
154 : {
155 : // Create the OpenMC mesh which will be tallied on.
156 722 : if (!_mesh_template_filename)
157 : {
158 : auto msh =
159 85 : dynamic_cast<const libMesh::ReplicatedMesh *>(_openmc_problem.getMooseMesh().getMeshPtr());
160 85 : if (!msh)
161 0 : mooseError("Internal error: The mesh is not a replicated mesh.");
162 :
163 : // Adaptivity and block restriction both require a map from the element subsets we want to
164 : // tally on to the full mesh.
165 85 : if (_use_dof_map)
166 : {
167 50 : _bin_to_element_mapping.clear();
168 :
169 : auto begin = _tally_blocks.size() > 0
170 50 : ? msh->active_subdomain_set_elements_begin(_tally_blocks)
171 70 : : msh->active_elements_begin();
172 50 : auto end = _tally_blocks.size() > 0 ? msh->active_subdomain_set_elements_end(_tally_blocks)
173 50 : : msh->active_elements_end();
174 182692 : for (const auto & old_elem : libMesh::as_range(begin, end))
175 182642 : _bin_to_element_mapping.push_back(old_elem->id());
176 :
177 : _bin_to_element_mapping.shrink_to_fit();
178 : }
179 :
180 : // When block restriction is active we need to create a copy of the mesh which only contains
181 : // elements in the desired blocks.
182 85 : if (_tally_blocks.size() > 0)
183 : {
184 40 : _libmesh_mesh_copy = std::make_unique<libMesh::ReplicatedMesh>(
185 20 : _openmc_problem.comm(), _openmc_problem.getMooseMesh().dimension());
186 :
187 40 : msh->create_submesh(*_libmesh_mesh_copy.get(),
188 40 : msh->active_subdomain_set_elements_begin(_tally_blocks),
189 40 : msh->active_subdomain_set_elements_end(_tally_blocks));
190 : _libmesh_mesh_copy->allow_find_neighbors(true);
191 : _libmesh_mesh_copy->allow_renumbering(false);
192 20 : _libmesh_mesh_copy->prepare_for_use();
193 :
194 20 : openmc::model::meshes.emplace_back(
195 40 : std::make_unique<openmc::LibMesh>(*_libmesh_mesh_copy.get(), _openmc_problem.scaling()));
196 : }
197 : else
198 65 : openmc::model::meshes.emplace_back(std::make_unique<openmc::LibMesh>(
199 65 : _openmc_problem.getMooseMesh().getMesh(), _openmc_problem.scaling()));
200 : }
201 : else
202 637 : openmc::model::meshes.emplace_back(
203 1274 : std::make_unique<openmc::LibMesh>(*_mesh_template_filename, _openmc_problem.scaling()));
204 :
205 722 : _mesh_index = openmc::model::meshes.size() - 1;
206 722 : _mesh_template = dynamic_cast<openmc::LibMesh *>(openmc::model::meshes[_mesh_index].get());
207 :
208 : // by setting the ID to -1, OpenMC will automatically detect the next available ID
209 722 : _mesh_template->set_id(-1);
210 722 : _mesh_template->output_ = false;
211 :
212 722 : _mesh_filter = dynamic_cast<openmc::MeshFilter *>(openmc::Filter::create("mesh"));
213 722 : _mesh_filter->set_mesh(_mesh_index);
214 722 : _mesh_filter->set_translation({_mesh_translation(0), _mesh_translation(1), _mesh_translation(2)});
215 :
216 : // Validate the mesh filters to make sure we can run a copy transfer to the [Mesh].
217 722 : checkMeshTemplateAndTranslations();
218 :
219 716 : return std::make_pair(openmc::model::tally_filters.size() - 1, _mesh_filter);
220 : }
221 :
222 : void
223 33 : MeshTally::resetTally()
224 : {
225 33 : TallyBase::resetTally();
226 :
227 : // Erase the OpenMC mesh.
228 33 : openmc::model::meshes.erase(openmc::model::meshes.begin() + _mesh_index);
229 33 : }
230 :
231 : Real
232 1263 : MeshTally::storeResultsInner(const std::vector<unsigned int> & var_numbers,
233 : unsigned int local_score,
234 : unsigned int global_score,
235 : std::vector<xt::xtensor<double, 1>> tally_vals,
236 : bool norm_by_src_rate)
237 : {
238 : Real total = 0.0;
239 :
240 1263 : unsigned int mesh_offset = _instance * _mesh_filter->n_bins();
241 2772 : for (unsigned int ext_bin = 0; ext_bin < _num_ext_filter_bins; ++ext_bin)
242 : {
243 1214821 : for (decltype(_mesh_filter->n_bins()) e = 0; e < _mesh_filter->n_bins(); ++e)
244 : {
245 1213312 : Real unnormalized_tally = tally_vals[local_score](ext_bin * _mesh_filter->n_bins() + e);
246 :
247 : // divide each tally by the volume that it corresponds to in MOOSE
248 : // because we will apply it as a volumetric tally (per unit volume).
249 : // Because we require that the mesh template has units of cm based on the
250 : // mesh constructors in OpenMC, we need to adjust the division
251 1213312 : Real volumetric_tally = unnormalized_tally;
252 1213312 : volumetric_tally *= norm_by_src_rate
253 1213312 : ? _openmc_problem.tallyMultiplier(global_score) /
254 949860 : _mesh_template->volume(e) * _openmc_problem.scaling() *
255 : _openmc_problem.scaling() * _openmc_problem.scaling()
256 : : 1.0;
257 1213312 : total += _ext_bins_to_skip[ext_bin] ? 0.0 : unnormalized_tally;
258 :
259 1213312 : auto var = var_numbers[_num_ext_filter_bins * local_score + ext_bin];
260 1213312 : auto elem_id = _use_dof_map ? _bin_to_element_mapping[e] : mesh_offset + e;
261 2426624 : fillElementalAuxVariable(var, {elem_id}, volumetric_tally);
262 : }
263 : }
264 :
265 1263 : return total;
266 : }
267 :
268 : void
269 722 : MeshTally::checkMeshTemplateAndTranslations()
270 : {
271 : // we can do some rudimentary checking on the mesh template by comparing the centroid
272 : // coordinates compared to centroids in the [Mesh] (because right now, we just doing a simple
273 : // copy transfer that necessitates the meshes to have the same elements in the same order). In
274 : // other words, you might have two meshes that represent the same geometry, the element ordering
275 : // could be different.
276 722 : unsigned int mesh_offset = _instance * _mesh_filter->n_bins();
277 401370 : for (int e = 0; e < _mesh_filter->n_bins(); ++e)
278 : {
279 400654 : auto elem_id = _use_dof_map ? _bin_to_element_mapping[e] : mesh_offset + e;
280 400654 : auto elem_ptr = _openmc_problem.getMooseMesh().queryElemPtr(elem_id);
281 :
282 : // if element is not on this part of the distributed mesh, skip it
283 400654 : if (!elem_ptr)
284 60654 : continue;
285 :
286 340000 : const auto pt = _mesh_template->centroid(e);
287 340000 : Point centroid_template = {pt[0], pt[1], pt[2]};
288 :
289 : // The translation applied in OpenMC isn't actually registered in the mesh itself;
290 : // it is always added on to the point, so we need to do the same here
291 : centroid_template += _mesh_translation;
292 :
293 : // because the mesh template and [Mesh] may be in different units, we need
294 : // to adjust the [Mesh] by the scaling factor before doing a comparison.
295 340000 : Point centroid_mesh = elem_ptr->vertex_average() * _openmc_problem.scaling();
296 :
297 : // if the centroids are the same except for a factor of 'scaling', then we can
298 : // guess that the mesh_template is probably not in units of centimeters
299 340000 : if (_openmc_problem.hasScaling())
300 : {
301 : // if scaling was applied correctly, then each calculation of 'scaling' here should equal 1.
302 : // Otherwise, if they're all the same, then 'scaling_x' is probably the factor by which the
303 : // mesh_template needs to be multiplied, so we can print a helpful error message
304 : bool incorrect_scaling = true;
305 147240 : for (unsigned int j = 0; j < OpenMCCellAverageProblem::DIMENSION; ++j)
306 : {
307 110430 : Real scaling = centroid_mesh(j) / centroid_template(j);
308 110430 : incorrect_scaling = incorrect_scaling && !MooseUtils::absoluteFuzzyEqual(scaling, 1.0);
309 : }
310 :
311 36810 : if (incorrect_scaling)
312 2 : paramError("mesh_template",
313 : "The centroids of the 'mesh_template' differ from the "
314 2 : "centroids of the [Mesh] by a factor of " +
315 2 : Moose::stringify(centroid_mesh(0) / centroid_template(0)) +
316 : ".\nDid you forget that the 'mesh_template' must be in "
317 : "the same units as the [Mesh]?");
318 : }
319 :
320 : // check if centroids are the same
321 : bool different_centroids = false;
322 1359992 : for (unsigned int j = 0; j < OpenMCCellAverageProblem::DIMENSION; ++j)
323 1019994 : different_centroids = different_centroids ||
324 : !MooseUtils::absoluteFuzzyEqual(centroid_mesh(j), centroid_template(j));
325 :
326 339998 : if (different_centroids)
327 4 : paramError(
328 : "mesh_template",
329 4 : "Centroid for element " + Moose::stringify(elem_id) +
330 8 : " in the [Mesh] (cm): " + _openmc_problem.printPoint(centroid_mesh) +
331 4 : "\ndoes not match centroid for element " + Moose::stringify(e) +
332 8 : " in the 'mesh_template' with instance " + Moose::stringify(_instance) +
333 8 : " (cm): " + _openmc_problem.printPoint(centroid_template) +
334 : "!\n\nThe copy transfer requires that the [Mesh] and 'mesh_template' be identical.");
335 : }
336 716 : }
337 : #endif
|
